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FOREWORD 


The  SPACE-A  program  described  in  this  report  is  a  modified  version 
of  the  trajectory  and  observation  generation  portion  of  the  Sequential  Position 
and  Covariance  Estimation  (SPACE)  program  originally  acquired  from  NASA. 
This  document  contains  a  detailed  analysis  of  the  equations  and  models  used 
by  the  program,  a  brief  description  of  the  routines  used  in  the  program,  as 
Well  as  instructions  for  its  operation.  The  principal  applications  of  SPACE-A 
are  the  prediction  of  space-vehicle  trajectories  and  the  generation  of  observa¬ 
tional  data  plus  other  trajectory  related  information. 

There  have  been  numerous  changes  made  in  the  original  program  in 
addition  to  its  adaptation  to  the  IBM  7030  and  its  subsequent  debugging.  The 
modified  SPACE-A  program  is  the  result  of  the  efforts  of  the  authors  and 
L.  E.  Wilkie.  S.  Schwartz  aided  in  the  debugging  and  rewriting  of  certain 
routines.  Special  thanks  and  appreciation  are  due  to  R.  K.  Squires,  D.  S. 
Woolston,  and  other  members  of  the  Special  Projects  Branch,  Theoretical 
Division  of  the  Goddard  Space  Flight  Center  from  whom  the  SPACE  program 
was  obtained. 

The  work  reported  in  this  document  was  performed  by  The  MITRE 
Corporation,  Bedford,  Massachusetts,  for  the  Directorate  of  Planning  and 
Technology,  Electronic  Systems  Division,  of  the  Air  Force  Systems  Command 
under  Contract  AF  19 (62 8) -5165. 


REVIEW  AND  APPROVAL 


Publication  of  this  technical  report  does  not  constitute  Air  Force  approval  of 
the  report's  findings  or  conclusions.  It  is  published  only  for  the  exchange  and 
stimulation  of  ideas. 

ANTHONY  P.  TRUNFIO 

Technical  Advisor,  Development  Engineering  Div. 

Directorate  of  Planning  and  Technology 
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ABSTRACT 

This  report  describes  the  SPACE-A  program  presently  available  for 
operational  use.  SPACE-A  is  the  trajectory  and  observation  generation 
portion  of  the  Sequential  Position  and  Covariance  Estimation  Program  (SPACE) 
which  is  currently  under  development.  The  document  contains  a  detailed 
analysis  of  the  equations  and  models  used  by  the  program,  a  brief  description 
of  routines  used  in  the  program,  as  well  as  instructions  for  its  operations. 

The  principal  applications  of  SPACE-A  are  the  prediction  of  the  space- 
vehicle  trajectories  and  the  generation  of  observational  data  plus  other 
trajectory  related  information. 
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SECTION  I 
OBJECTIVES 


INTRODUCTION 

The  SPACE  trajectory  detenninatlon  program  was  originally  devel¬ 
oped  for  the  Special  Projects  Branchy  Theoretical  Division,  Goddard 
Space  Flight  Center  by  the  Sperry  Rand  Systems  Group.  It  was  designed 
as  a  comprehensive  trajectory  determination  and  tracking  program  for 
both  orbital  and  deep  space  flights.  The  original  program  consisted 
of  three  major  modes  of  operation: 

(1)  SPACE-A,  designed  for  trajectory  and  observation  generation 
without  statistical  processing. 

(2)  SPACE-Bl,  designed  for  statistical  estimation  of  the  six 
state  variables  describing  the  position  and  velocity  at 
prescribed  points  of  the  trajectory, 

(3)  SPACE-B2,  designed  for  statistical  estimation  of  the  six 
state  variables  describing  position  and  velocity  plus  up 
to  20  additional  states  selected  from  a  group  of  vari¬ 
ables  which  permit  estimation  of  dynamic  biases  (para¬ 
meters  affecting  orbital  motion)  or  observational  biases 
(affecting  measurements). 

The  version  of  the  SPACE  program  received  by  MITRE  contained 
numerous  errors  in  programming,  SPACE-Bl  and  SPACE-B2  had  never  been 
debugged  and  many  of  the  routines  of  SPACE-A  were  unworkable.  In  ad¬ 
dition,  the  programming  had  to  be  made  compatible  with  the  IBM  7030, 
Therefore,  a  large  effort  was  put  into  checking  the  theory,  debugging 
the  program,  and  revising  or  rewriting  certain  routines.  Many  check¬ 
out  runs  were  necessary  during  and  after  the  debugging  of  the  program. 

Since  the  SPACE-B2  mode  includes  the  capabilities  of  the  SPACE-Bl 
mode,  the  SPACE-Bl  mode  was  not  debugged.  The  original  SPACE-B2 
(which  was  not  operational  when  received  by  MITRE)  employed  two  dif¬ 
ferent  statistical  estimation  techniques: 

(1)  a  minimum  variance  sequential  filtering  technique,  the  so- 
called  ’*Kalman  filter*',  and 


1 


(2)  a  non-recurslve  batch  processing  technique* 

Because  a  number  of  trajectory  estimation  programs  at  MITRE 
already  employ  the  batch  processing  technique,  the  second  option  of  the 
SPACE-B2  mode  was  left  in  a  non-opcrational  status*  The  first  option 
of  SPACE-B2^  which  incorporates  the  Kalman  filter  in  the  trajectory 
estimation  procedure,  is  presently  being  debugged  and  checked  out* 

The  theory,  programming,  and  results  of  its  operation  will  be  pub¬ 
lished  in  a  succeeding  document.  The  theoretical  development  and 
much  of  the  programming  of  the  trajectory  generation  portion  of 
SPACE-B2  is  identical  with  that  used  in  SPACE-A. 

This  document  deals  with  the  modified  SPACE-A  trajectory  genera¬ 
tion  program  written  in  Fortran  IV  and  is  presently  operational  on 
the  IBM  7030  at  MITRE.  This  report  is  based  on  the  original  documents^ 
published  by  Sperry  Rand  for  NASA. 

Section  I  contains  a  brief  discussion  of  the  objectives  and  capa¬ 
bilities  of  SPACE-A.  Section  II  discusses  the  theoretical  background, 
as  well  as  the  equations  and  methods  employed  by  the  program  Section 
III  consists  of  a  user's  guide  for  operating  the  program.  Section  FV 
gives  the  program  structure  and  a  brief  description  of  the  function 
of  each  subroutine. 

The  larger  part  of  the  report  is  contained  in  Section  II  since 
understanding  or  following  the  programming  is  really  a  matter  of  fol¬ 
lowing  the  appropriate  equations.  Although  parts  of  Section  TI  are 
original,  the  contents  of  pertinent  sections  of  the  SPACE  Analytic 
Manual^ have  been  freely  used  or  modified  to  fit  the  description 
of  the  present  program.  Sections  III  and  IV  closely  follow  the  format 
of  the  original  User's  Manual and  Programmer's  Manual  however, 

since  a  number  of  changes  in  programming  have  been  made,  only  Section 
III  and  Section  IV  of  this  report  should  be  used  in  operating  the  modi¬ 
fied  SPACE-A  trajectory  generation  program. 
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CAPABILITIES 


The  racxiified  SPACE-A  trajectory  generation  may  be  run  in  any  of 
the  following  modes: 

(1)  Trajectory  generation  -  the  computation  of  position  and 
velocity  of  the  vehicle  at  regular  time  points  along  a 
trajectory. 

(2)  Observation  generation  -  the  computation  of  certain 
ground-based  observations  of  the  vehicle  at  regular 
time  points. 

(3)  Simulated  data  -  the  computation  and  writing  onto 
tape  of  the  observations  in  a  format  suitable  for 
processing  data  by  SPACE-B2. 

(4)  Visibility  computation  -  the  time  of  initial  view 
of  the  vehicle  by  a  given  station,  the  total  time 
in  sight  by  the  station,  and  the  observations 
while  the  vehicle  is  in  sight. 

The  required  input  consists  of  the  initial  date  and  initial  time 
of  a  particular  run  and  the  initial  position  and  velocity  of  the  vehi¬ 
cle.  When  the  effect  of  drag  is  to  be  included  which  is  usually  the 
case,  the  effective  area  and  mass  of  the  vehicle  and  the  magnitude  of 
solar  flux  must  be  specified.  When  the  effect  of  direct  solar  radia¬ 
tion  is  included,  an  effective  area  pertaining  to  solar  radiation  must 
be  specified.  For  observation  calculations,  the  specification  of  the 
station  location  on  the  surface  of  the  Earth  is  necessary.  The  de¬ 
tails  of  the  units  and  format  of  these  required  inputs  along  with 
other  optional  input  features  are  discussed  in  Section  III.  The  types 
of  output  data,  which  consist  primarily  of  trajectory  and/or  observa¬ 
tion  information,  is  also  given  in  Section  III. 

The  philosophy  behind  the  modified  SPACE-A  trajectory  generation 
program  is  that  its  primary  use  would  be  for  Earth  orbital  trajec¬ 
tories.  Toward  this  objective,  its  dynamic  model  includes  the  effects 
of  the  primary  terra  and  oblateness  perturbations  of  the  Earth’s  grav¬ 
ity,  the  solar,  lunar,  and  planetary  gravitational  perturbations,  the 
effect  of  direct  solar  radiation  and  drag,  as  well  as  such  dynamic  ef¬ 
fects  as  precession,  nutation,  and  the  daily  rotation  of  the  Earth. 
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The  model  used  for  drag  is  a  combination  of  the  U*S.  Standard 
Atmosphere  1962  and  the  Harris-Priester  model  for  the  upper  atmosphere, 
which  depends  on  the  magnitude  of  solar  flux*  The  model  for  Earth 
gravitational  oblateness  may  include  up  to  nine  zonal  harmonics,  four 
sectorial  harmonics,  and  twelve  additional  tesseral  harmonics.  The 
types  of  ground  station  observations  that  may  be  specified  are  given 
in  Section  II  (see  OBSERVATIONS). 

Although  the  modified  SPACE-A  program  is  intended  primarily  for 
Earth  orbital  trajectories,  it  may  be  used  for  other  types  of  missions. 
However,  a  number  of  options  that  were  included  in  the  original  SPACE-A 
trajectory  generation  program  are  not  presently  operational.  These  op¬ 
tions  were  a  capability  for  modeling  the  perturbations  due  to  thrust, 
a  model  for  lunar  libration,  simple  models  for  lunar  gravitational  ob¬ 
lateness  perturbations,  as  well  as  models  for  the  drag  of  Mars  and 
Venus,  and  a  capability  for  computing  on-board  observations  from  a 
vehicle.  Although  these  options  are  not  presently  available  in  the 
modified  SPACE-A  program,  with  some  effort  they  could  also  be  included. 
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SECTION  II 
THEORY 


GENERAL  DESCRIPTION  OF  TRAJECTORY  COMPUTATION 

Orbit  prediction  or  trajectory  computation  is  the  process  of  cal¬ 
culating  the  position  and  velocity  of  a  vehicle  at  any  time  subsequent 
to  some  initial  time,  provided  the  initial  position  and  velocity  of 
the  spacecraft  are  given. 

To  accomplish  this  prediction,  one  uses  the  laws  of  celestial 
mechanics  embodied  in  the  differential  equations  of  motion.  The 
forcing  functions  for  the  equations  of  motion  are  obtained  from  a  dy¬ 
namic  model  which  accounts  for  the  accelerations  on  the  vehicle.  A 
reference  frame  is  erected  to  express  the  components  of  the  various 
vector  quantities,  and  the  equations  of  motion  are  numerically  inte¬ 
grated  subject  to  the  given  initial  conditions.  Once  the  vehicle  tra¬ 
jectory  is  determined,  the  program  can  then  generate  observations. 

The  coordinate  system  used  in  this  program  is  based  upon  the 
position  of  the  mean  Earth’s  equator  and  equinox  at  0^  1  January  of 

the  year  subsequent  to  the  initial  input  time  of  the  program.  Coordi¬ 
nate  directions  of  this  frame  are  inertial  with  respect  to  the  fixed 
stars;  the  center  of  origin  of  the  system,  however,  may  be  transferred 
from  one  central  body  to  another  so  that  the  spacecraft  motion  is 
specified  relative  to  a  point  mass  which  itself  has  a  proper  motion. 
This  reference  frame  is  called  the  Base  Date  System  or  simply  the  in¬ 
ertial  system. 

Observations  made  from  the  Earth  are  necessarily  in  a  system 
different  from  the  Base  Date  System.  To  accomplish  transfer  between 
various  coordinate  systems,  transformations  are  provided  (see  COORDI¬ 
NATE  SYSTEMS  AND  TRANSFORMATIONS).  These  transformations  include  the 
dynamical  effects  of  precession,  nutation,  and  daily  rotation  of  the 
Earth. 
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All  accelerations  acting  on  the  vehicle  are  specified  In  the  Base 
Date  System*  The  gravitational  attractions  of  bodies  In  the  solar  sys¬ 
tem  are  functions  only  of  position  with  respect  to  the  vehicle;  conse- 
quently^  the  program  employs  an  ephemerls  giving  planetary  coordinates 
relative  to  the  Sun^  and  lunar  coordinates  relative  to  the  Earth,  all 
In  a  Base  Date  System*  A  Base  Date  System  Is  specified  for  overlapping 
two-year  blocks  of  data,  the  date  corresponding  to  the  middle  of  the 
two-year  file*  Specifying  an  Initial  time  causes  the  program  to  choose 
an  ephemerls  file  having  as  Its  Base  Date  the  beginning  of  the  year 
following  the  Initial  Input  time  of  the  program.  In  this  way,  at  least 
one  full  year  of  ephemerls  Information  Is  available  before  a  change  of 
reference  system  Is  necessary* 

Another  acceleration  specified  In  the  Base  Date  System  without 
transformation  Is  that  arising  from  solar  radiation  pressure*  Since 
this  acceleration  Is  a  function  of  relative  position  between  the  Sun 
and  the  vehicle.  Its  direction  Is  given  In  the  proper  frame  by  using 
information  from  the  ephemerls* 

Other  accelerations,  such  as  perturbations  due  to  Earth  oblate¬ 
ness  effects  must  be  transformed  through  nutation  and  precession  to  the 
proper  frame*  All  positions,  velocities  and  accelerations  are  computed 
In  canonical  units  (i*e*.  Earth  radii  (ER) ,  Earth  radii  per  hour 
(ER/hr),  Earth  radii  per  hour  per  hour  (ER/hr^),  respectively)  and  ap¬ 
propriate  constants  are  used  to  transform  to  other  units  of  measure. 

Vehicle  motion  Is  always  computed  relative  to  some  reference  body; 
a  planet;  the  Moon;  the  Sun*  Consequently,  the  equations  of  motion 
contain  a  term  which  accounts  for  the  acceleration  of  the  reference 
body  on  the  vehicle.  The  remaining  accelerations  are  usually,  but  not 
always,  much  smaller  than  this  primary  acceleration  and  are  therefore 
called  perturbations*  In  most  cases,  they  can  be  regarded  as  giving 
rise  to  small  disturbances  in  the  orbit  determined  by  the  reference 
body  and  the  Initial  conditions* 
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Reference  bodies  are  changed  during  a  trajectory  calculation  when 
the  vehicle  leaves  the  **reglon  of  Influence”  associated  with  a  partic¬ 
ular  body.  Regions  of  Influence  are  computed  for  a  body  with  respect 
to  the  object  of  which  it  is  a  satellite.  Hence,  each  planet  has  a 
region  of  influence  defined  relative  to  the  Sun,  and  the  Moon  h*is  a 
similar  region  defined  relative  to  the  Earth,  In  transferring  into  or 
out  of  such  a  region,  velocity  as  well  as  position  with  respect  to  the 
new  reference  body  must  be  calculated. 

Since  no  analytic  solution  exists  for  the  equations  of  motion, 
numerical  methods  are  employed  to  compute  the  components  of  position 
and  velocity.  In  the  program,  a  choice  may  be  made  between  using 
straightforward  integration  and  using  Encke's  method.  The  former 
technique,  called  Cowell's  method,  is  conceptually  simple,  but  suffers 
from  precision  and  machine  running  time  problems,  Encke's  method, 
although  somewhat  more  complicated,  gives  dividends  in  both  precision 
and  machine  efficiency.  In  this  procedure,  the  Keplerlan  orbit 
arising  from  the  reference  body  central  force  field  is  taken  as  a 
nominal  trajectory.  The  perturbation  accelerations  are  Integrated 
and  the  resulting  position  and  velocity  Increments  are  added  to  the 
Keplerlan  solutions.  Naturally,  Encke's  method  is  most  effective  when 
the  perturbations  are  small. 

The  present  program  has  the  capability  to  compute  a  number  of 
ground-based  observations •  Corrections  are  provided  for  the  refraction 
of  an  electromagnetic  signal  by  the  troposphere  and  by  the  ionosphere. 
Adjustments  are  computed  for  errors  in  elevation  angle,  range,  and 
radial  velocity;  other  angular  corrections  are  calculated  from  the 
adjustment  in  elevation  angle. 

Description  of  the  basic  equations  used  by  the  program  for  the 
dynamic  model  now  follows.  It  includes  derivations  of  the  accelera¬ 
tions  due  to  the  planets.  Earth  oblateness,  drag,  and  direct  solar 
radiation  pressure.  It  also  discusses  the  coordinate  systems,  numeri¬ 
cal  integration  methods,  and  observation  calculations  Including  the 
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refraction  model  used  by  the  program# 

DYNAMIC  MODEL 

The  equations  of  motion  for  a  space  vehicle  are  second-order 
differential  equations  which  describe  the  accelerations  arising  from 
the  forces  acting  on  the  vehicle.  These  forces  can  be  classified  as 
follows ; 

(1)  Gravitational  acceleration  of  the  reference  body  -  primary. 

(2)  Gravitational  perturbations  due  to  other  bodies  (e,g., 
planets) . 

(3)  Gravitational  perturbations  due  to  reference  body 
oblateness . 

(4)  Perturbations  due  to  drag. 

(5)  Perturbations  due  to  direct  solar  radiation  pressure 
on  the  space  vehicle. 

For  orbit  determination  the  reference  body  primary  gravitational 
force  almost  always  predominates  the  perturbation  forces  given  above. 
The  only  exception  occurs  during  reentry  of  a  vehicle  into  the  atmos¬ 
phere,  when  drag  forces  may  be  larger  than  the  primary  force  of  gra¬ 
vity.  Another  force  which  may  exceed  the  primary  force  of  gravity  is 
that  of  vehicle  thrust.  Since  this  program  is  designed  primarily  for 
orbital  computations,  this  force  is  not  discussed  here  although  it 
can  be  readily  included  if  desired. 

The  simplest  gravitational  force  field  is  that  due  to  a  single 
point  mass  or  equivalently  that  due  to  a  homogeneous  ponderable  body 
which  is  perfectly  spherical.  In  this  case  the  equations  of  motion 
of  a  vehicle  with  respect  to  the  ponderable  reference  body  are: 


where 

P 

G 

M 

m 


G(M+m)  =  GM,  since  m  <<  M 
is  the  universal  gravitational  constant, 
is  the  mass  of  the  reference  body, 
is  the  mass  of  vehicle. 


(1) 
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R  is  the  vector  position  of  the  vehicle  with  respect  to 
the  reference  body  center. 

With  initial  conditions  Rq  and  Rq  Equation  (1)  defines  a  **two- 
body"  or  Keplerian  orbit,  which  arises  from  the  primary  reference 
body  central  force  field  and  which  may  be  expressed  in  closed  form  in 
terms  of  its  true  or  eccentric  anomaly. 

The  above  "two-body**  dynamic  model  may  be  extended  to  Include 
per turbational  accelerations  acting  on  the  vehicle,  in  which  case 
Equation  (1)  becomes: 


(2) 


where 


Pj  is  the  summation  of  all  the  accelerations  arising  from 
planetary, lunar  and  solar  attraction, 

P2  is  an  acceleration  arising  from  the  oblateness  of  the 
reference  body, 

P3  is  an  acceleration  arising  from  drag  as  derived  from 
an  appropriate  model  to  be  described  later,  and 

P4  is  an  acceleration  due  to  direct  solar  radiation  upon 
the  vehicle  neglecting  reflected  sunlight  from  the 
reference  bodies  or  planets. 

The  perturbational  accelerations  Pj,  P2,  P3,  and  P4  are  obtained 
from  specific  dynamical  models  whose  details  are  described  below. 

PERTURBATIONS  DUE  TO  PLANETARY  ATTRACTIONS 

The  general  expression  for  the  perturbational  acceleration,  Pj, 
of  a  space  vehicle  due  to  the  gravitational  influence  of  the  Sun,  Moon, 
and  other  planets  (excluding  the  reference  body)  is  given  by: 


j 


(3) 


where 


yj  -  GMj 

G  is  the  gravitational  constant. 


is  the  mass  of  the  jth  body. 
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Ryj  is  the  poBition  of  the  vehicle  with  respect  to  the 
jth  body,  (Figure  1) 

Rrj  is  the  position  of  the  reference  body  with  respect 
to  the  jth  bodyi  (Figure  1) 

AR  is  the  position  of  the  vehicle  with  respect  to  the 
reference  body  (Figure  1) • 


jth  body 


Figure  1.  Vectors  for  Planetary  Attractions 


If  the  two  terms  of  the  bracketed  expression  in  Equation  (3)  are 
very  nearly  equals  a  loss  of  accuracy  due  to  round-off  errors  by 
machine  computations  will  occur.  Equation  (3)  can  be  written  in  a 
more  convenient  form  due  to  Battin^^^  thereby  eliminating  this  pro¬ 
blem.  Using  Equation  (3) ^  and  Equation  (4)  below  and  rearranging  terms 
one  obtains  Equation  (5). 

AR  -  Ryj  -  Rrj  (4) 


This  latter  form  is  used  to  achieve  more  computational  accuracy.  The 
actual  programmed  equations^  which  can  readily  be  shown  equivalent  by 
substitution  are: 

Pi  "  I  t^rj  (f(U))  -  AR] 

where 

f(U)  -  Ut3  -m  (3-HU).] 

1  +  (1+U)  ^ 


(5) 


(6) 


(7) 


u 


2 (Rrj  +  AR)  •  ^ 


(8) 
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The  forces  Introduced  by  solar,  lunar  and  planetary  attractions 
in  orbit  prediction  about  the  Earth  are  usually  smaller  by  a  factor 
of  10“^  than  that  of  the  primary  attraction  of  the  Earth’s  gravita¬ 
tional  pull*  Planetary  perturbations  (including  solar  and  lunar  per¬ 
turbations)  do  have  a  significant  effect  on  a  vehicle  trajectory  over 
long  time  periods  or  for  deep  space  probes* 

PERTURBATIONS  DUE  TO  EARTH  OBLATENESS 

The  fact  that  the  Earth  is  not  a  perfect  sphere  of  uniform  den¬ 
sity  gives  rise  to  perturbational  accelerations  due  to  Earth  oblate¬ 
ness.  Formulas  are  derived  below  which  the  program  uses  to  include 
these  perturbations  in  the  prediction  of  near-Earth  trajectories.  In 
a  coordinate  system  attached  rigidly  to  the  Earth,  the  Earth  oblate¬ 
ness  perturbation  may  be  treated  as  the  acceleration  of  a  conservative 
force  derivable  from  a  potential.  Thus: 

P2  -  -  VU  (9) 

The  potential  function  is  given  by  Equation  (12)  and  the  final  form  of 
the  acceleration  ^2  given  by  Equations  (26)-(29)  and  Equation 
(32)  below. 

The  Geopotential  Function 

The  geopotential  function  can  be  derived  from  the  basic  property 
that  the  potential,  U,  satisfies  Laplace’ s  Equation^ 


V^u  -  0* 


(10) 


In  spherical  coordinates.  Equation  (10)  becomes 


V^u 


W  -L  r 


1  +  _L  f  1  m.'l 

94)^  9X  ^cos  ({>  9X' 


0  (11) 


Although  the  solution  of  Equation  (11)  is  not  elementary,  the  equa¬ 
tion  may  be  solved  by  applying  the  method  of  separation  of  variables 
and  using  Legendre  polynomials^^’  jhg  solution  of  Equation 
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(12) 


(11)  is  given  by: 


n-o 


in«o 


where 


U  is  the  Earth's  gravitational  parameter » 

Re  is  the  Earth's  mean  equatorial  radius » 

are  spherical  coordinates  (see  Figure  2), 

Pn  is  the  associated  Legendre  function  (spherical  harmonics) 
of  the  first  kind  of  degree  n  and  order  m,  and 

Cn^m  and 

Sn^m  are  numerical  coefficients. 

The  Legendre  polynomials,  defined  by 


and  the  associated  Legendre  function 


Pn(x)  -  (1  -  x2)7  JSL  Pn(x). 


dx^ 


The  spherical  coordinates  r.  A,  (J)  are  defined  with  respect  to 
the  fundamental  planes  determined  by  the  true  equator,  and  the  Green¬ 
wich  meridian  (see  Figure  2). 


z 


Greenwich  meridian 


vehicle 


^  y 


true  equator 


X 


Figure  2.  Spherical  Coordinates 
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The  fundamental  term  in  the  expression  for  U  is  given  by 
m  ■  n  ■  0.  The  terms  in  which  m  ■  0  are  called  zonal  harmonics. 
Inspection  of  Equation  (12)  shows  that  these  terms  vary  only  with 
latitude,  and  hence  reflect  deviations  of  the  Earth's  potential  from 
a  sphere  of  uniform  density  that  are  symmetric  around  the  spin  axis 
(e.g.,  a  pear  shaped  Earth  potential  can  be  modeled  with  zonal  harmon¬ 
ics).  Evaluating  the  2,0  term  for  U  we  have 

-  ^  (sin  (^)  C2,o 

and  since 

P2(x)  -  I  (3x2  _  1) 

Equation  (13)  becomes 

"  r  2  ^2,0. 

The  zeros  of  Equation  (13)  are  at  *  35?25. 


Figure  3.  Zonal  Harmonic  for  n  =  2,  m  =  0 

The  sectorial  harmonics  arise  when  m  ■  n.  The  zeros  of  the 
sectorial  terms  lie  along  lines  of  longitude.  The  remaining  combi¬ 
nations  are  tesseral  or  square  in  that  these  terms  vanish  both  along 
a  number  (m  -  n)  of  parallels  of  latitude  and  a  number  (2m)  of  meri¬ 
dians  of  longitude.  Figure  4  provides  an  illustration  of  zonal, 


(13) 
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V 


sectorial,  and  tesseral  harmonic  variations  for  a  sample  set 


Zonal 

Harmonics 


Sectorla,! 

Harmonics 


Tesseral 

Harmonics 


m 


n 


-  4 

-  0 


n  •  3 
m  ■  3 


n 


-  8 


Figure  4,  Variations 


In  summary,  Equation  (12)  describes  the  Earth's  gravitational 
potential  at  any  point  in  space;  the  negative  gradient  of  this  poten¬ 
tial  gives  the  corresponding  acceleration.  Note  should  be  made,  how¬ 
ever,  that  the  acceleration  obtained  is  in  the  Greenwich  coordinate 
system,  and  hence  must  be  transformed  into  the  inertial  system  (the 
Base  Date  System)  to  be  employed  in  trajectory  computation. 

Computation  of  Acceleratlon^Dus-ttr^garth  Oblateness 

The  components  of  gravitational  acceleration  are  now  expressed 
as  the  sum  of  terms  which  have  the  same  general  form  for  any  m,  n 
combination.  First,  it  is  necessary  to  develop  expansions  for  cos  mX 
and  sin  mX  in  terms  of  cos  X  and  sin  X,  From  DeMoivre's  theorem. 


cos  mX  +  i  sin  mX  ■  (cos  X  +  i  sin  X)^ 


Now,  expanding  by  the  binomial  theorem 


For  m  an  even  number 


14 


cos 


a 

2 


mX  +  i  sin  mX  -  I  ( (-l)'^(co8  X)"'"‘^^(sin  X) 


.m-2k. 


2  k 


k«0 


'a_i 

2 


k,  ,.m-2k-l,  .  ,v2k+l  > 
(cos  X)  (sin  X)  j  » 


k*o 


and  for  m  an  odd  number 


m-1 

2 


cos  mA  +  1  sin  mX  -  J]  (^j^)  (-l)^(cos  X)^"'^^(sin  X) 


vm-2k , 


?> 


2 


k-0 


m  ,xk,  __  ,,in-2k-l,__ 


^  ^  ^2k+lk"^) 


k-o 


(sin  X)' 


Equating  real  and  imaginary  parts  above,  we  have 


M 


cos 


k“o 

M' 


mX  -  I  (-1)'^  (  ^)  (cos  X)"’"'^*^  (sin  X) 


m-2k 


.2k 


sin  m\  ~  I  (-1)  (2^+1^*^°*  **  ""^(sin  X)‘ 


k-o 


rrp 


where  (“)  - 


m! 


if  m  is  even 


M 


if  m  is  odd  i 


"  1  if  ®  i®  even 


and  M’ 


m^l 


if  m  is  odd  a 


From  the  definition  of  spherical  coordinates: 
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sin  4) 


cos  4)  cos  i  ■  —  / 
r 


cos 


4)  sin  i 


(16) 


where  and  z  are  the  coordinates  of  the  vehicle  in  the  Green¬ 

wich  coordinate  system. 

Using  Equation  (16)  in  Equations  (14)  and  (15),  we  have: 


cos  mX  ■ 


1. 


M 


(r  cos  (}>)“  k-o 


r  ,  -.kcmi,  ,.m-2k,  ,.2k 

I  (-1) 


(17 


M' 


sin  mX  ■ 


m-2k-l .  V 2k+l 

(y) 


(r  cos  (J))  k-o 

Substitution  of  Equations  (17)  and  (18)  into  Equation  (12)  yields: 

Il-fj  I  y) 


n-0 


m-o  (r  cos  (J>) 


where 


M 


G(x,  y)  -  Cn,m  I  (-D^l^j^.)  (x) 
k-0 


m-2k,  v2k 

(y) 


(]S) 


(12') 


M' 


+  s„.„  I 

k-0 


in-2k-l,  ,.2k+l 

(y) 


From  the  definition  of  the  associated  Legendre  function,  we  have 


m 


P^T) 


X) 

2"  n! 


-  t2  J 


dx 


m+n 


[(t2  -  1)"]. 


m 


Setting  T  ■  sin  ^  and  dividing  Equation  (19)  by  cos  0, 


(39) 


Pn  (sin  ») 

cos^  4) 
n 


2‘*  n!  dx 


^  [(T^  -  1)"1 
m+ti 


Setting  U  -  I  Ujn^ni  using  Equation  (20)  we  have 

n*0  m^O 


(20) 
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Utn.n 


»“  r 


' 

JL_. 

n 

1 

dtiH-n 

^  r  ^ 

.r®  2"  nl. 

,  tn+n 

dT 

(t^  -  1)^  0(x,  y) 


-UR, 


m+n+l  , 
r  2  n! 


I  m+n 


di 


trrfn 


(t2  -  I)""  G(x,  y) 


The  gravitational  acceleration  is  computed  by  taking  the  negative 
gradient  of  the  potential  function.  The  general  derivative  in  the 
gradient  will  be  taken  with  respect  to  where  5  takes  on  the 

values  of  x,  y^  and  z. 

Defining 


au 


n>m 


and  carrying  out  the  indicated  operations  using 


H  r 


and 


ii .  ^  [i]  i  i  ris.  _  ^] 

35  35  LrJ  r  [H  r2j 


yields : 


M  Rg" 


^n,n 


2" 


dT^ 


_  G(x,_  y)  [iz  _  z5l  (t2  -  1)'^ 

r  135  r2j  ,  m+n+l 


(t2  -  1) 


dx 
n  aG(x 


dT 


nH-n 


y)1 

H  J 


where 


-  Cn,m  I  [(-1)''  (2"  )f("»“2k)  y^'^ 


35 


k-0 


H 


,  m-2k  2k-l  Byi-, 

+  2k  X  y 


M  * 

r  r/,\k^  m  •)(/  ,x  in“2k— 2  2k+l  3x 
+  Sn.in  I  [(-1)  (2vliJ{(ni-2k-l)  X  y  — 


k-0 


35 


in-2k-l  2k  ^ 
35' 


+  (2k+l)  x^  "*'  "  y*--'  l^}]. 


(21) 


(22) 
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that 


where 


It  can  be  ahown  by  applying  the  binomial  theorem  and  Induction 

4  -  D"  ■"  I  ^'(-1)“  (”.) 

k»0 

[n  -  4] 


and 


■j]  Indlcatsa  the  integral  part  of  Hence, 


jm+n 


[221] 


(  2  _  .  y  (“■)  /  \n-2k-m 

^  m+n  ^  ^  IjcJ  (n-2k-m)! 

dx  k-o 


iflLiili  (t2  -  d"  -  I  (-l)'^  ["]  (^)n-2k-in-l 

n,+r.+1  L  u;  fn-2k-m-n! 


m-n-1 , 

I  2  J 


,  m+n+1 
dT  k-o 


k'  (n-2k-m-l)! 


(23) 


(2A 


where  t  -  sin  ^  .  Therefore,  Equation  (21)  may  be  rewritten  as 


[2?i] 


^n,m 


”  ^  ^9  ft  "T*  ,  ,.k  .ni  (2n-2k)  !  cz.^n-2k-mj 

-n  ,  mfn+1  L‘  ^  ^  ^k^  (n-2k-m) !  4-’  • 

2  n!  r  k-O 


{Sslspi  5  G(x.  y)  -  - 


3^ 


1  I  (S)  a;r|» 


k-0 


(25) 


By  separately  considering  the  case  for  n-m  odd  and  for  n-m  even, 
it  may  be  shown  that  Equation  (25)  reduces  to  the  following  expres¬ 


sion: 


[“] 

u  ^  2  J 


.C  V  f,  ,.k  _ _  rz^n-2k-m-l,  ^ 

An.m  -  -  I  k!  (nV) !  (n-2k-m) !  4^ 


2  r 


k-0 


{((2n-2k.H)  ^  c(,,  y)  -  f 
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where 


(27) 


2x  _  rl  if  C  -  X  ^  _  rl  if  C  -  y  iz  rl  if  C  - 

3C  "  '0  if  £  X*  3£  "  'O  if  £  y*  3£  "  'O  if  £  ?< 


and 


M 


G(x»  y)  - 


^Ti^XR  I  (“1) 
k-0 


jn!  ^  xra-2k,  ^  2k 

(x)  (y) 


2k I  (m-2k)! 


M' 


,m-2k-l,  ,2k+l 


.  o  V  /  _ SJ _  /  xm-zk-j 

+  Sn,ni  I  (-1)  (2k+l)!  (in-2k-l)  ! 

k-0 


(2R) 


M 


m/2  if  m  is  even 


m-l 


M’  - 


Y  “  1  if  m  is  even 


If  m  Is  odd 


tn-1 


if  m  is  odd 


M 


I  I(-l) 


ml 


k-0 


2k  I  (m-2k)l 


{ (m-2k) (x) 


m-2k-l,  v2k  3x 

‘y'  — 


+  2k(x)'"-^Ny)=''‘-l  |2)1 

9£ 


M' 


Snjm  I  [ (*1) 


ml 


k-0 


(2k+l)!  (m-2k-l)! 


//  i\/  vin-2k-2.  v2k4*l  9x  .  /ni^-i\/  Nni-2k-l,  .2k  9y i , 

X  {(m-2k-l)(x)  (y)  rr  +  (2k+l)  (x)  (y)  ^|] 


(29) 


It  is  important  to  note  that  when  z  -  0,  a  special  procedure 
must  be  used  since  the  term  arises  in  Equation  (26)  for  n-m 

even;  and  when  n-m  is  odd. 

When  n-m  is  even.  Equation  (26)  reduces  to  one  term,  when  z  =  0. 

n^m 

'(n4m4-l)£  G  ^ 


_  ~  ^  -  f(  1)  ^  (irHn)  1  M  Kn4in+]J 


(30) 


Similarly,  when  n-m  is  odd.  Equation  (26)  may  be  reduced 


19 


An,m 


0  for  C  “  K|  y 


Re 

nH-n+l 

2  r 


(^)  (n-[»])l/ 


£ 


(31) 


In  summary 9  the  program  computes  the  acceleration  of  a  vehicle 
due  to  the  Earth  In  the  Greenwich  coordinate  system  by  the  formula: 


P2  -  I  I  t  +  A^.m  J  +  k)  (32) 

n-0  m"0 

where  jf  k  are  unit  vectors  In  the  Greenwich  coordinate  system. 

In  Equation  (32)  the  range  of  indices  (n,ra)  are  restricted  by 
SPACE-A  to  the  following  limits; 

n  m 

2  0,  1,  2 

3  0,  1,  2,  3 

4  0,  1,  2,  3,  4 

5  0,  1,  2.  3 

6  0,  1,  2 

7  0,1 

8  0 

9  0 

10  0 

The  fundamental  term  is  given  by  the  (0,  0)  combination  and 
would  be  the  only  term  present  if  the  Earth  were  of  uniform  density 
and  a  perfect  sphere*  Since  the  center  of  the  coordinate  system  is 
taken  to  be  at  the  center  of  mass,  it  can  be  shown  that  the  (1,  0) 
and  (1,  1)  combinations  vanish^ ^  . 

F 

In  Equation  (32),  A^  ^  is  given  by  (26)  when  z  ^  0,  and  when 
z  =  0,  by  Equation  (30)  or  Equation  (31)  depending  on  whether  n-m 
is  even  or  odd. 
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PERTURBATION  DUE  TO  DRAG 


Accurate  simulation  of  an  artificial  satellite  or  other  space 
vehicle  trajectories  requires  consideration  of  vehicle  deceleration 
resulting  from  atmospheric  drag.  A  number  of  planets  (e*g.^  Earth, 
Venus,  Mars  and  Jupiter)  have  sufficiently  dense  atmospheres  to  re¬ 
tard  the  motion  of  a  vehicle  within  their  atmospheres;  however,  this 
discussion  confines  itself  to  the  atmospheric  model  of  the  Earth. 

An  analysis  of  drag  must  take  into  account  the  particular  mis¬ 
sion  of  the  vehicle,  e.g.,  low  eccentricity,  orbit,  reentry,  or  fly¬ 
by,  since  vehicle  mission  determines  what  portion  of  the  atmosphere 
it  is  necessary  to  include  in  the  model. 

The  following  discussion  describes  the  general  equations  used 
for  drag  computation,  some  of  the  problems  involved  in  simulating  the 
Earth’s  atmosphere,  and  the  effects  of  certain  simplifying  assumptions 

Drag  Equations 

The  program  uses  two  different  formulas  for  computing  the  vehi- 
cle  deceleration,  P3,  due  to  drag.  For  a  relatively  dense  atmos¬ 
phere  where  the  assumption  of  continuum  flow  is  valid  the  following 
well  known  equation  is  used: 

-►  _1  CpS  ^ 

P3  -  -  Cl  P  — ]  Va  Va 

where 

p  is  the  atmospheric  density  at  the  vehicle, 

CD  is  the  drag  coefficient  of  the  vehicle  (dimensionless), 

S  is  the  effective  surface  area  of  the  vehicle, 

m  is  the  mass  of  the  vehicle, 

Va  is  the  vector  of  the  velocity  of  the  vehicle  with 
respect  to  the  local  atmosphere, 

Va  is  the  magnitude  of  Va,  and 

Cl  is  a  constant  used  to  convert  the  above  expression 
into  the  units  used  by  the  program. 

Equation  (33)  is  used  to  compute  the  acceleration  of  drag  in  the 


(33) 
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lower  atmosphere  from  0  kilometers  up  to  120  kilometers,  although  the 
model  for  the  lower  atmosphere  from  which  the  values  of  p  and 
are  obtained  may  be  extended  up  to  210  kilometers  with  some  loss  of 
accuracy.  The  basic  model  used  in  the  lower  atmosphere  is  the  U.  S , 
Standard  Atmosphere  1962^^^  and  its  details  are  discussed  belov;. 

As  the  atmosphere  becomes  more  diffuse,  the  mean  free  path  length 
(average  distance  between  impacts  of  air  molecules)  increases.  At 
110  kilometers,  mean  free  path  length  is  roughly  one  meter  and  at  130 
kilometers  may  become  as  large  as  ten  meters^ When  the  mean  free 
path  length  exceeds  the  dimensions  of  the  vehicle,  the  assumption  of 
continuum  air  flow  is  no  longer  applicable.  In  such  a  diffuse  atmos¬ 
phere  where  all  collisions  are  essentially  two-body  collisions,  the 
air  flow  is  referred  to  as  free  molecular  flow. 

Ketchum^^^  has  derived  the  following  formula  for  the  magnitude 
of  drag  deceleration  in  free  molecular  flow: 

IP3I  -  i  (1  +  TltP  Cav  f]  Va 

where 

R  the  radius  of  the  vehicle, 

X  the  mean  free  path, 

Cav  the  average  velocity  of  particles  in  the  medium. 

Ketchum  is  uncertain  as  to  the  validity  of  the  (1  +  2R/X)  term 
in  Equation  (34).  In  the  high  atmospheric  region  the  assumption  that 
X  >>  R  is  usually  justified,  except  perhaps  below  140  kilometers. 
Therefore,  the  program  actually  uses  Equation  (35)  in  computation  of 
the  drag  in  the  high  atmosphere. 

P3  ■  -  C2  P  ^av  *]  Va 

where, 

S  is  the  effective  surface  area  of  the  vehicle, 
m  is  the  mass  of  the  vehicle, 
p  is  the  atmospheric  density  at  the  vehicle, 

Cav  is  the  average  velocity  of  particles  in  the  medium. 


(34) 


(35) 
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Vg  is  the  vector  of  the  velocity  of  the  vehicle  with 
respect  to  the  surrounding  atmosphere |  and 

Cp  is  a  constant  used  to  convert  the  expression  into 
the  units  used  by  the  program. 

Equation  (35)  is  used  to  compute  the  deceleration  of  drag  in  the 
upper  atmosphere  from  100  kilometers  up  to  1^050  kilometers,  although 
the  model  for  the  upper  atmosphere  may  be  extended  up  to  2,000  kilo¬ 
meters.  The  basic  model  used  for  the  upper  atmosphere  is  that  due  to 

r9i 

Harris  &  Priester*”  .  The  values  of  o  and  Cav  used  in  Equation 
(35)  arc  derived  from  this  model,  the  details  of  which  are  also  rlis- 
cussed  below. 

Notice  in  both  Equations  (33)  and  (35)  that  the  direction  of  the 
drag  force  is  in  a  direction  opposite  to  the  velocity  with  respect  to 
air.  In  addition,  both  formulations  assume  zero  lift  and  assume  that 
the  angle  of  attack  of  the  vehicle  is  zero,  i.e,,  that  the  vehicle 
velocity  relative  to  the  air  mass  is  in  line  with  the  vehicle  longi¬ 
tudinal  axis. 

It  is  readily  seen  that  the  formula  for  drag  in  the  upper  atmos¬ 
phere,  Equation  (35),  differs  from  the  equation  of  drag  in  the  lower 
atmosphere,  Equation  (33),  Furthermore,  in  the  region  between  100 
kilometers  and  120  kilometers  there  are  disagreements  between  the  two 
models.  For  example,  the  value  of  density  predicted  by  the  low  atmos¬ 
pheric  model  (U,  S,  Standard  Atmosphere  1962)  at  ]20  kllorcters  is 
34%  higher  than  that  predicted  by  the  high  atnosplieric  model  (-’arris- 
Pricstcr)  at  the  same  height. 

The  present  program  achieves  a  compromise  between  these  t^'^’o 
models  by  treating  the  region  from  100  to  120  kilometers  as  n  trr-^s-*- 
tion  region  in  which  a  weighted  average  is  taken  between  the  drag 
values  computed  by  the  two  methods  and  gradually  slides  the  weight 
from,  unity  for  free  molecular  flow  and  zero  for  continuum  flow  at  120 
kilometers  to  unity  for  continuum  flow  and  zero  for  free  molecular 
flow  at  100  !;ilometers. 
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Because  of  the  uncertainties  in  the  atmospheric  models  and  be¬ 
cause  of  the  approximations  made  in  the  analysis  the  computation  of 
drag  deceleration  is  probably  accurate  to  +  5%  in  the  lower  atmosphere 
and  is  less  accurate  in  the  upper  atmosphere. 

Variables  Used  in  Drag  Computation 

Vehicle  Mass 

In  the  most  general  case^  the  vehicle  mass  in  the  drag  equations 
should  be  considered  as  variable  with  time.  In  the  orbiting  case  or 
the  fly-by  case,  a  step  change  in  mass  representing  the  separation  of 
a  landing  craft  is  conceivable,  A  long-term  steady-state  mass  flow 
rate,  however,  would  probably  be  small. 

For  the  reentry  case,  if  the  reentry  vehicle  is  of  the  heat-sinK 
type,  the  mass  would  be  constant.  For  an  ablative  nose  cone  ^i.e., 
one  which  loses  mass  when  moving  at  high  speeds  due  to  friction)  the 
mass  flow  rate  is  a  function  of  the  drag.  For  ballistic  missil<^  ap¬ 
plications,  this  mass  change  is  usually  ignored.  In  any  event,  such 
changes  in  mass  represent  a  small  error  in  the  location  of  the  impact 
point.  Therefore,  the  program  treats  the  mass  as  an  input  constant. 

Surface  Area 

The  effective  surface  area  term  S  in  the  drag  equation  is  not 
simply  the  cross-sectional  area  of  the  vehicle.  The  vehicle,  in 
passing  through  the  air,  produces  a  shock  wave  which  skirts  the  mis¬ 
sile  thereby  placing  the  effective  cross-sectional  area  at  a  point 
somewhat  close  to  the  nose.  Since  the  shock  wave  changes  with  air 
speedy  so  does  the  effective  cross-sectional  area.  In  practice,  S 
is  made  constant  and  any  variation  with  speed  is  included  in  the  co¬ 
efficient  of  drag. 

Velocity  with  Respect  to  Air 

The  velocity  of  the  vehicle  is  available  in  an  inertial  coordi¬ 
nate  system  (Base  Date  System), 

The  vehicle  velocity  with  respect  to  the  moving  air  mass,  Va, 
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in  the  same  coordinate  systemi  is  obtained  by  subtracting  the  velocity 
of  the  air  mass  from  the  vehicle  velocity#  A  good  first  approximation 
to  the  velocity  of  the  air  mass  is  obtained  by  assuming  the  air  mass 
to  be  rigidly  attached  to  the  rotating  planet# 

From  these  considerations  we  have; 

Va  -  R  -  n'  X  R 

where 

is  the  velocity  with  respect  to  surrounding  atmosphere ^ 

“► 

R  is  the  position  vector  of  the  vehicle  in  the  inertial 

frame. 

#  ^ 

R  is  the  velocity  vector  of  the  vehicle  in  the  inertial 
frame 9  and 

n  is  the  vector  of  angular  rotation  of  Earth  expressed 
in  the  inertial  frame# 

A  better  approximation  could  be  obtained  by  including  the  effects 
of  wind  velocity#  The  purely  local  effects  have  to  be  neglected »  but 
the  long-term  horizontal  effects  are  known  as  a  function  both  of  posi¬ 
tion  on  the  Earth's  surface  and  of  altitude#  The  effects  of  the  wind 
velocity's  direction  (independent  of  altitude  but  dependent  on  lati¬ 
tude  and  longitude)  and  magnitude  (strongly  dependent  on  altitude » 
less  strongly  on  latitude^  and  least  on  longitude)  would  have  to  be 
included#  The  error  made  by  neglecting  Earth  winds  is  about  1,500 
feet  at  impact  for  a  typical  ICBM  mission#  It  should  be  noted  that 
winds  are  of  importance  only  in  the  Earth's  lower  atmosphere,  mainly 
for  the  reentry  case#  In  the  present  program  the  effect  of  winds  is 
neglected# 

Drag  Coefficient 

The  drag  coefficient,  Cd,  is  sometimes  considered  to  be  a  con- 

CdS 

stant#  Very  often  the  ballistic  coefficient,  B  -  — —  is  used  in 

"  m 

analysis  in  place  of  the  drag  coefficient,  mass,  and  effective  sur¬ 
face  area#  A  much  more  accurate  representation  for  Cd  is  obtained 
by  considering  it  to  be  a  function  of  Mach  number#  Mach  number  is 


(36) 
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defined  to  be: 


where 9 

M  Is  the  Mach  number ^ 

|V|^|  Is  the  speed  with  respect  to  the  surrounding  alr^  and 
c  Is  the  speed  of  sound  In  the  surrounding  alr« 

The  speed  of  sound  Is  a  function  of  altitude  obtained  from  the 
low  atmospheric  model  (see  below).  It  should  be  noted  that  as  altitude 
Increases,  the  atmosphere  becomes  rarlfled  to  the  point  that  the  speed 
of  sound  loses  Its  physical  significance. 

In  the  program  Cd  Is  tabulated  for  about  AO  different  Mach 
numbers.  These  numbers  are  denser  for  speeds  below  Mach  2  than  those 
above  and  are  very  dense  In  the  region  around  Mach  1,  For  Intermedi¬ 
ate  values  of  Mach  number,  linear  Interpolation  Is  used. 

Inadequate  knowledge  of  the  drag  coefficient  Is  one  of  the  major 
sources  of  Inaccuracy  In  the  simulation  of  drag.  Since  drag  coeffi¬ 
cient  Is  a  function  of  Mach  number,  drag  coefficient  data  has  been 
obtained  by  wind  tunnel  measurements  made  at  a  range  of  Mach  numbers. 

Air  Density 

In  the  region  below  120  kilometers  the  air  density,  p»  Is  a 
function  of  altitude  obtained  from  the  low  atmospheric  model  and  Is 
computed  from  a  stored  table  (see  below) , 

In  the  high  atmosphere  p  Is  obtained  from  the  upper  atmospheric 
model  of  Harrls-Prlester  and  Is  considered  to  be  a  function  of  alti¬ 
tude,  local  solar  time,  solar  flux,  and  latitude.  Tables  are  provided 
In  the  program  for  Its  computation  (see  below). 

Mean  Particle  Velocity 

The  mean  particle  velocity,  Cav»  i®  of  Importance  only  when 
the  vehicle  Is  In  the  upper  atmosphere  where  the  assumption  of  free 
molecular  flow  Is  valid*  From  Equations  I*3«A-1  and  X«2,6-l  of 
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(38) 


Reference  [7],  is  given  by 


where 

T  is  the  absolute  temperature  of  particles  (*K), 
m  Is  the  mean  molecular  weight  of  medium  (gm) , 
k  is  a  constant  of  proportionality  (.145)# 

The  values  of  T  and  m  are  given  in  the  Harris-Priester  model  of 
the  upper  atmosphere.  Cav  can,  therefore,  be  obtained  directly  from 
stored  tables. 

Lower  Atmospheric  Model 

Drag  in  the  lower  atmosphere  (below  120  kilometers)  can  be  large 
and  a  vehicle  entering  this  region  will  usually  be  slowed  down  suffi¬ 
ciently  to  be  quickly  captured  by  the  Earth,  Thus,  the  lower  atmos¬ 
phere  is  primarily  of  concern  in  the  reentry  case. 

Data  for  an  average  model  have  been  well  established  for  the 
lower  atmosphere.  There  are  five  sources  for  these  data:  U.  S, 
Standard  Atmosphere  1962:  COSPAR  International  Reference  Atmosphere 
(CIRA),  1961;  COESA  Table  for  Tropical  Latitudes,  1962;  ARDC  Model 
Atmosphere  1956,  1959,  Table  I  shows  the  density  deviation  (in 
percent),  as  a  function  of  altitude,  of  each  of  the  others  from  the 
U.  S.  Standard  Atmosphere  values.  From  the  table,  it  is  evident  that, 
except  for  the  COESA  tables,  there  is  good  agreement  between  the  vari¬ 
ous  tables  at  low  altitudes.  Note  that  the  U.  S.  Standard  Atmosphere 
and  CIRA  tables  are  in  excellent  agreement  all  the  way  to  120  km, 
(400,000  feet). 

The  lower  atmosphere  is  characterized  by  seasonal,  diurnal,  and 
latitude  variations;  however,  none  of  these  is  sufficiently  well 
documented.  The  only  effect  of  omitting  them  is  that  the  impact  point 
of  a  re-entering  body  would  be  slightly  different.  It  was  estimated 
in  1958  that  the  standard  deviation  for  a  heat-sink  type  nose  cone 
used  in  the  ICBM  application  is  only  about  0,5  nm. 
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Table  I 


Comparison  of  Sources  of  Density  Data 


Al' 

L i tude 

U.  S.  Standard 
Atmosphere 
Density 

Percent 

Deviatio- 

r ror  I'f  : 

- 1 

♦ 

: .  rc  ’■( 

Values 
(Reference) 
slugs/ f t^ 

ARDC 

1956 

ARDC 

1959 

PTP  A 

corsA 

1962 

kn. 

ft. 

Lr  1  I\A 

1961 

0 

0 

2.38"3 

0 

0 

0.55 

-4.77 

3.0 

10,000 

1.76"3 

0 

0 

-0.91 

-5.32 

5.5 

18,000 

1.36“3 

0.04 

0 

1.85 

-1.67 

ll'.l 

33,000 

7.97“^* 

0.05 

0 

1.68 

1.92 

lA .  b 

48,000 

4.00"^* 

0.09 

0 

2.36 

15.5 

20.4 

67,000 

l.bl”*" 

3.28 

0.16 

0.48 

6.80 

29.0 

95,000 

4.20"5 

0.59 

-2.36 

0.10 

0.46 

32.5 

110,000 

2.07"5 

-3.13 

-3.13 

0.68 

2.53 

4b. b 

160,000 

2.32"^ 

4.77 

4.77 

0.77 

8.9J 

b7.1 

220,000 

2.50"'^ 

15.0 

15.5 

1.30 

8.10 

31.4 

300,000 

4.62-^ 

31.2 

-10.8 

0.11 

-  - 

1  121.9 

400,000 

3.62“^ 1 

81.5 

-35.0 

1.17 

"  J 

Three  tables  with  32  values  each  are  stored  in  the  program  for 
the  lower  atmospheric  model.  The  first  contains  32  values  of  alti¬ 
tudes  in  kilometers  I  the  second  contains  32  values  of  the  logarithm 
(base  10)  of  the  density,  p,  in  gm/km^,  and  the  third  contains  32 
values  of  the  speed  of  sound,  c,  in  km/sec.  The  n^^  entry  in  the 
altitude  table  corresponds  to  the  n^^  entry  in  the  table  for  logjoP 
and  c.  Intermediate  values  are  obtained  by  linear  interpolation. 
These  three  tables  are  obtained  from  the  U,  S,  Standard  Atmos¬ 
phere  1962,  in  which  densities  and  speed  of  sound  at  altitudes  below 
90  km  are  listed.  Above  90  km  the  speed  of  sound  was  calculated  from 
temperature  and  mean  molecular  weight  data  which  are  directly  avail¬ 
able. 
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The  value  of  altitude  (in  Earth  radii)  above  an  ellipsoidal 
Earth  is  obtained  using  a  formula  found  in  Baker^^^^. 

h  -  r  -  1  +  f  sin2  +  -y  (“  -  "I)  sin^  2^  (39) 

where 

h  is  the  height  of  the  vehicle  (ER) , 

r  is  the  distance  of  the  vehicle  from  the  geocenter  (ER) , 
is  the  geocentric  latitude  of  the  vehicle, 

f  is  the  flattening  constant  of  the  Earth  -  Vq q  ' V  • 


For  programming  purposes  Equation  (39)  is  put  into  the  more  convenient 


form: 


h  -  R  -  1  + 


~  2^U98.3 


298.3  R2 


(40) 


where , 

h  is  the  altitude  of  the  vehicle  (ER) , 

X|y,z  are  the  position  coordinates  of  the  vehicle  in  the 
Greenwich  coordinate  system  (ER), 


2  is  the  distance  of  the  vehicle  from 
the  geocenter  (ER) • 

The  position  coordinates  of  the  vehicle  are  available  from  the  pro¬ 
gram  and  h  is  multiplied  by  6378.165  to  convert  its  units  into 
kilometers. 

After  h  is  obtained,  p  (p  ■  10^^®^®^)  and  c  are  obtained 
from  the  low  atmosphere  tables.  Then  the  air  velocity,  Va  and  the 
drag  coefficient,  Cp,  are  computed  as  described  above.  Finally, 

P3  is  computed  according  to  Equation  (33). 

Upper  Atmospheric  Model 

Models  of  the  Earth* s  upper  atmosphere  (above  100  kilometers) 
must  take  into  account  solar  activity.  There  is  evidence  that  solar 
activity  occurs  cyclically  at  periods  of  27  days,  6  months,  1  year 
and  11  years. 


R  ~  /I 


+  +  z 
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Theoretical  models  do  not  exist  for  the  27-dayi  6-month|  and  1- 
year  cycles*  Diurnal  varlatlonSp  If  anyi  of  the  models  for  these 
cycles  are  not  known*  Investigation  of  the  11-year  cycle  (corres¬ 
ponding  to  the  sunspot  period)  in  solar  flux  led  to  the  Harris- 

1  9  111 

Priester  model  of  the  upper  atmosphere.  This  model*"  ’  has 
diurnal  and  solar  flux  variations. 

r9i 

The  Harrls-Prlester  model  lists  the  density,  absolute  tem¬ 

perature,  and  mean  molecular  weight  of  the  atmosphere  as  a  function 
of  altitude,  solar  flux  and  local  solar  time*  The  mean  particle 
velocity  can  be  found  by  use  of  Equation  (38).  Note  that  at  the 
North  and  South  Poles  local  solar  time  Is  undefined. 

The  upper  atmosphere  has  a  delaying  effect  on  solar  radiation* 

It  takes  several  hours  for  the  Sun's  heat  to  pass  through  the  atmos¬ 
phere  and  reach  the  Earth's  surface*  The  Harrls-Prlester  model  Is 
based  on  densities  computed  at  the  Earth's  equator*  Intuitively,  It 
Is  expected  that  It  will  take  longer  for  the  solar  flux  to  reach  the 
poles  as  opposed  to  the  equator*  Therefore,  It  Is  considered  that 
there  Is  an  effective  variation  of  solar  flux  with  latitude.  This 
variation  Is  Implemented  In  the  program  by  applying  the  Harrls-Prlester 
model  at  the  equator  and  a  stored  table  of  "twilight”  densities  at 
the  poles.  The  cosine  of  the  latitude  of  the  vehicle  Is  used  as  a 
weighting  factor  to  Interpolate  between  the  two  sets  of  data* 

The  Harrls-Prlester  upper  atmospheric  model  has  been  Incorporated 
In  the  program  by  means  of  a  table  look-up  procedure*.  Two  tables  are 
used  In  the  program* 

The  first  table  lists  the  logarithm  to  the  base  10  of  the  density 
In  gm/km^  times  the  mean  particle  velocity  In  km/sec  (log^oP  ^av) 
the  equator  as  functions  of  altitude  (km),  solar  flux  (watts/m^  -  Hz 
at  a  wavelength  of  10*7  cm)  and  local  solar  time  (hrs).  These  values 
have  been  tabulated  for  16  values  of  altitude,  4  values  of  solar  flux 
and  13  values  of  local  solar  time*  Hence,  this  table  has  832  entries* 
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The  second  table  lists  logjo  PCav  In  the  same  units  at  **twi- 
llght”  for  the  modeling  of  the  polar  region  as  a  function  of  altitude 
and  solar  flux*  This  table  has  16  entries  for  altitude »  A  entries 
for  solar  fluX|  or  6A  entries. 

For  Intermediate  values  of  the  Input  variables  (altitude,  solar 
flux,  local  solar  time)  a  linear  Interpolation  Is  used  to  obtain  the 
output,  logio  P^av*  This  method  gives  fairly  accurate  results  since 
oCav  Is  nearly  exponential.  The  value  of  solar  flux  Is  determined  by 
Input  data.  The  value  of  altitude  Is  computed  according  to  Equation 
(AO),  The  value  of  local  solar  time  Is  computed  from  the  x  and  y 
coordinates  of  the  vehicle  and  the  Sun  In  the  Inertial  (Base  Date) 
coordinate  system  (see  Figure  5), 


Figure  5.  View  of  x-y  Plane  from  Above 

Thus  we  have: 

local  solar  time  ■{(02  -  0i  +  tt)  4.  720®}  (mod  360®) 


where 

ei 

02 


tan"^  (“) 
tan"^  [^] 


where  ^  01  ^ 
where  ^  i 
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Xy,  yv  are  the  vehicle  coordinates  in  the  inertial  system, 

^8»  ys  the  Sun's  coordinates  in  the  inertial  system* 

Local  solar  time  given  above  is  then  divided  by  15  to  convert  its  units 
to  hours. 

After  a  value  of  logjo  pCav  bas  been  obtained  from  both  the 
equatorial  table  and  the  '^twilight**  polar  table  an  interpolation  is 
done  to  obtain  the  final  value  of  pCav  using  the  cosine  of  latitude 
as  a  weighting  function*  Hence: 

L  «  Lp  +  cos  ^  (Le  “  Lp) 


2  +  y^  + 
pCav  “  10 

where 

L  is  the  final  value  of  logjo  PCav» 

Lp  is  the  value  of  log^o  POav  obtained 
from  the  twilight  table, 

Le  is  the  value  of  logjo  PCav  obtained 
from  the  equatorial  table, 

x,y,z  are  the  coordinates  of  the  vehicle  in  the  inertial 
(Base  Date)  system  and  ^  is  the  latitude  of  the 
vehicle. 

Once  pCav  has  been  found,  the  vehicle  area  and  mass  (S  and  m, 
standard  program  inputs)  as  well  as  air  velocity  (see  Equation 

(36)),  are  obtained  and  the  value  of  drag  deceleration  is  finally 
computed  according  to  Equation  (35). 

PERTURBATION  DUE  TO  DIRECT  SOLAR  RADIATION 

Solar  radiation  exerts  a  pressure  on  the  intercepting  surface  of 
a  vehicle.  Orbiting  planetary  vehicles,  having  a  large  area  to  mass 
ratio  are  subject  to  noticeable  perturbations  due  to  solar  radiation 
pressure.  In  fact,  for  orbits  above  500  miles  the  solar  radiation 


(43) 

(44) 
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perturbation  is  usually  more  significant  than  that  due  to  drag^ 

The  vehicle  acceleration  due  to  solar  radiation  pressure  depends  on 
the  area  to  mass  ratio  of  the  vehicle  as  well  as  the  intensity  of  the 
Sun’s  incident  power  at  the  vehicle  and  the  fraction  of  solar  illu¬ 
mination  on  the  vehicle*  The  illumination  factor  is  considered  in 
three  distinct  regions  in  the  following  analysis:  full  sunlight^ 
penumbral  illumination^  and  no  illumination  (i.e*^  umbral  region). 

In  computing  the  solar  radiation  perturbation,  P4,  this  analy¬ 
sis  neglects  the  dispursive  effects  of  planetary  atmospheres  which 
complicate  the  geometry  of  the  umbra  and  penumbra.  The  analysis  also 
neglects  the  effects  of  reflected  sunlight  from  the  reference  body 
or  any  other  planet. 


Acceleration  Pue  to  Solar  Radiation  PreasurB 


An  expression  for  the  acceleration  due  to  solar  radiation  pres- 

[121 

sure  given  in  Wolverton  is: 


where. 


p  is  the  illumination  factor, 

q  is  the  reflectivity  coefficient, 

A  is  the  area  of  vehicle  pertaining  to  solar 
radiation  pressure, 

m  is  the  mass  of  vehicle, 

Lo*3.86  X  10^^  watts,  the  total  power  output 
of  the  Sun  (+  3%) , 

c  is  the  speed  of  light, 

Rsv  is  the  vector  from  the  Sun  to  the  vehicle. 

The  SPACE  program  assumes  a  reflectivity  coefficient  equal  to 
unity.  This  assumption  may  be  altered  by  changing  the  area.  A,  by 
an  appropriate  factor.  Simplifying  Equation  (45)  and  putting  the 
above  variables  into  units  used  by  the  program,  one  obtains: 


(45) 


(46) 
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where > 

p  is  the  illumination  factor,  0  <  p  ^  1, 

A  is  the  area  pertaining  to  radiation  pressure  (ft^), 

m  is  the  mass  pertaining  to  radiation  pressure  (lb-masses), 

Rsv  is  the  vector  the  Sun  to  the  vehicle  (ER) , 

Cp  -1,04819  X  10^  ®  constant  used  to  convert 


the  expression  into  units  used  by  the  program  and  to  ac- 
LO 

count  for  the  factor  of  Equation  (45). 

Equation  (46)  is  the  basic  equation  used  by  the  program  for  com¬ 
puting  the  solar  radiation  perturbation  acceleration. 

The  illumination  factor,  p,  is  obtained  from  the  relative  geo¬ 
metry  of  the  vehicle  with  respect  to  the  Sun,  the  Hocn  and  the  planet*" 
p  -  1  in  full  sunlight, 
p  -  0  in  the  umbral  region, 

0  <  p  <  1  in  the  penumbral  region. 

It  is  possible  that  a  vehicle  may  lie  within  the  penumbral  region 
of  two  bodies  at  the  same  time,  e,g,,  the  Earth  and  the  Moorv  In 
this  case,  only  the  penumbral  illumination  factor  due  to  the  reference 
body  is  computed*  Errors  introduced  by  this  assumption  are  extremely 
small:  first,  because  for  most  vehicles  of  interest  the  solar  radia¬ 
tion  perturbation  is  itself  small  (usually  less  than  10  ^  times  the 
acceleration  of  the  reference  body  except  for  low  density  balloons) ; 
second,  because  only  a  short  time  is  spent  in  the  penumbral  region  of 
the  reference  body  by  an  orbiting  vehicle;  and  third,  because  the 
incidence  of  simultaneous  penumbral  obscuration  is  rare. 

Since  the  geometry  used  in  calculating  the  illumination  factor 
is  the  same  as  that  used  for  eclipse  information,  the  portion  of  the 
program  which  computes  the  solar  radiation  pressure  perturbation  also 
computes  the  times  at  which  the  vehicle  enters  or  leaves  the  umbra  or 
penumbra  of  the  celestial  bodies.  The  geometry  used  in  calculating 
the  vehicle  illumination  factor,  p,  is  described  below. 
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Vehicle  Illumination  Factor 


Figure  6  illustrates  the  geometry  of  vehicle  Illumination,  neg¬ 
lecting  the  effects  of  atmospheric  refraction.  By  similar  triangles 
it  is  seen  that  the  height  of  the  umbral  cone  (hy)  and  penumbral 
cone  (hp)  are  respectively  given  by: 


where. 


Rsp  is  the  distance  from  the  center  of  the  Sun  to  the 
center  of  the  reference  body, 

Rs  is  the  radius  of  the  Sun,  and 

Rp  is  the  radius  of  the  reference  body. 

Next,  criteria  are  developed  to  see  in  which  region  the  vehicle 
lies,  i.e.,  full  sunlight,  umbra,  or  penumbra.  First,  consider 
Figures  7  and  8  and  the  definitions  and  relations  below. 

Rgp  is  the  vector  from  the  center  of  the  Sun 
to  the  center  of  the  reference  body, 

R  is  the  vector  from  the  center  of  the 
reference  body  to  the  vehicle. 


hu 


Rsp 


(A7) 


(A8) 


(A9) 


Q  -  -  hp 

^  Ksp 


From  Figure  7  we  get  the  following: 


(50) 


cos  A  -  /l  - 


COS  a 


(R  -  P)  »  ^8P 

1^  -  p|  Rsp 


(51) 

(52) 


35 


36 


Figure  6.  Vehicle  Illumination 


POSITION  OF  VEHICLE 


APEX  OF  UMBRA 


Figure  7.  Umbral  Region  Geometry 


POSITION  OF 


Figure  8.  Penumbrol  Region  Geometry 
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Tccsi  Figur%3  8  we  get\ 


cos 


(53; 


CO.  6  -  »,  -  %  ■.  (54; 

1^  -  Q!  Rap 

If  the  scalar  product  then  the  vehicle 

lies  on  the  side  of  the  reference  body  away  from  the  Sun,  In  this 
case,  if  cos  a  >  0  and  if  ot  ^  A,  the  vehicle  or  satellite  lies 
in  the  unbia  and  p  is  set  equal  to  zero.  Another  test  is: 

if  I  cos  a  I  21  1  cos  a] 

then  p  *  0| 


i.e.p  the  vehicle  lies  in  the  umbra. 

If  the  vehicle  does  not  lie  in  the  umbral  region  one  can  test  to 
see  if  it  lies  in  the  penumbra,  Thus^  If  cos  6  >  0  and  6  ^  B  the 
vehicle  lies  in  the  penumbra,  or  equivalently i 

if  I  cos  6|  21  I  cos  b|  (5^) 

then  0  <  p  <  1 

j^i.e,,  the  vehicle  lies  in  the  penumbra). 

If  the  vehicle  does  not  lie  in  either  the  umbral  or  penumbral 
region,  then  the  reference  body  does  not  obscure  the  Sun's  rays,  A 
check  can  then  be  made  to  see  whether  any  other  celestial  body  blocks 
or  partially  blocks  the  solar  radiation;  and  if  not,  the  illumination 
factor  is  set  equal  to  one,  p  ■  1, 


Penumbral  Illumination  Factor 

If  the  tests  of  Equations  (55)  and  (56)  above  indicate  that  the 
vehicle  lies  in  the  penumbral  region,  then  the  illumination  factor, 
p,  is  computed  according  to  the  formulat 


P 


Be* 


(57; 
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where  Aex  Is  the  enguler  aree  subtended  by  the  exposed  portion  of 
the  solar  cap  at  the  vehicle  position^  and 

6s  Is  the  total  angular  area  of  the  solar  cap  at  the 
vehicle  position. 

The  general  approach  used  by  the  SPACE  program  for  computing  the 
penumbral  Illumination  factor  Is  somewhat  more  detailed  than  that 
given  In  most  references.  It  consists  In  projecting  the  solar  disc 
(or  cap)  and  the  cap  of  the  obscuring  planet  (or  moon)  on  to  a  great 
imaginary  sphere  whose  center  Is  at  the  vehicle  position.  The  rela¬ 
tive  angular  areas  of  the  caps  are  computed  and  the  Illumination  fac¬ 
tor  Is  given  according  to  Equation  (57).  While  the  theoretical  de¬ 
velopment  Is  somewhat  Involved ^  It  leads  to  simple  closed  form  alge¬ 
braic  expressions  convenient  for  use  on  a  computer. 

Firsts  consider  Figure  9  which  shows  a  sphere  of  radius  re¬ 

presenting  the  Sun  or  a  planet »  and  a  point  P  at  a  distance  t  ^ 
from  the  center  of  the  sphere.  The  apparent  angular  area  of  the 
sphere  as  seen  from  point  P  Is  equivalent  to  the  angular  area  of 
the  sphere's  cap  projected  onto  a  projection  sphere  of  radius  a. 
Therefore  the  following  Is  truei 


angular  area  6  ■  2Tr  [1  -  cos 


■  2if 


.  Aa  +  2R) 

L  i  +  R  . 


where  a  Is  the  radius  of  the  projection  sphere,  and 
H  Is  the  height  of  the  projected  cap. 

The  total  angular  area  of  a  planet  (6p)  or  of  the  Sun  (65)  at 
the  position  of  a  vehicle  can  be  obtained  by  substituting  the  appro¬ 
priate  values  of  I  and  R  Into  Equations  (58)  and  (59)  and  simpli¬ 
fying.  Thus,  we  have: 


(58) 

(59) 


0p  ■  2ir 


/h(h  +  2Rp) 


1  - 


h  +  Rp 


-  2Tt 


I’ 


(60) 
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POSITION  OF  VEHICLE  AT 


Figure  9.  Geometry  Used  to  Measure  the  Angular  Area 
of  the  Solar  or  Planetary  Cop 
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I 


(61) 


where 

6p  Is  the  angular  area  of  the  planet  (or  moon)  at  the  vehicle 
position^ 

h  Is  the  height  of  the  vehicle  above  the  planet's  surface^ 

Rp  Is  the  radius  of  the  planet ^ 

6$  Is  the  angular  area  of  the  Sun  at  the  vehicle  position^ 

Rsv  Is  the  distance  from  the  center  of  the  Sun  to  the  vehicle^ 

Rs  Is  the  radius  of  the  Sun^ 

H  Is  the  height  of  the  solar  cap^ 

H'  Is  the  height  of  the  planetary  cap^ 

a  Is  the  radius  of  the  Imaginary  sphere  of  projection. 

Equations  (60)  and  (61)  establish  the  relative  size  of  the  solar  and 
planetary  caps  as  seen  from  the  vehicle. 

Next 9  It  Is  desired  to  determine  the  angular  area  of  the  exposed 
solar  cap  Agx*  *^0  do  thls^  consider  Figure  10.  In  this  figure  the 
vehicle  Is  positioned  at  the  center  of  a  great  Imaginary  projection 
sphere  of  arbitrary  radius^  a.  The  relative  positions  of  the  Sun 
and  the  obscuring  body  or  planet  are  shown  as  well  as  the  projections 
of  the  solar  and  planetary  caps  onto  the  great  sphere.  The  equator 
of  the  great  sphere  Is  constructed  to  be  coplanar  with  the  vehicle^ 
the  center  of  the  Sun,  and  the  center  of  the  planet.  A  great  circle 
Is  also  constructed  perpendicular  to  the  equator  and  passing  through 
the  Intersection  of  the  caps;  this  circle  will  be  used  to  define  one 
of  the  limits  of  Integration  In  the  computation  of  the  angular  area 
of  the  two  lunes  which  are  formed  by  the  great  circle. 

For  the  calculations  to  follow,  three  coordinate  systems  are  de¬ 
fined  and  Illustrated  In  Figures  10  and  11.  The  coordinate  systems 
employed  are; 


A1 


SOLAR 

TOTAL 


POSITION  OF  VEHICLE 


and  Planetary  Caps 


y$  y'  t  y” 


Figure  11  .  Definition  of  Coordinate  System 


(1)  A  syatam  (x,y,z)  for  the  solar  cap. 

(2)  A  system  (x'.y'.z')  for  the  planetary  cap. 

(3)  A  system  (x",y",z")  for  the  great  circle. 

All  three  coordinate  systems  are  orthogonal  and  have.  In  common,  the 
same  y  axis.  The  relations  between  the  coordinates  are  given  In 
Equations  (62)  and  (63) t 


X  " 

cos  0c 

sin  6c 

0 

X 

z' 

- 

sin  6c 

cos  0c 

0 

Z 

_y’j 

0 

0 

1 

y 

"x"' 

sin  6q 

cos  6g 

o" 

X 

z" 

- 

-  cos  6c 

sin  6q 

0 

z 

y" 

0 

0 

1 

y 

where 

6c  Is  the  angle  from  the  center  of  the  solar  cap  (z-axls)  to 
the  center  of  the  planetary  cap  (z'-axls), 

6g  Is  the  angle  from  the  center  of  the  solar  cap  (z-axls)  to 
the  great  circle  (x"-axls). 

By  Inspection  of  Figure  10  one  finds  that: 


(62) 


(63) 


cos  dc 


R«v 


(64) 


where 

0  <  6c  i  Y 

Next  find  the  points  of  intersection  of  the  planetary  and  solar 
caps  X  •  Xp*  y  ■  yp^  and  z  •  zp#  Referring  to  Figures  9  and  10 
one  sees  that  the  equation  of  the  small  circle  of  the  solar  caps 
(having  a  cap  height  of  H)  ist 
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z  -  a  -  H  (65) 

x2  +  y2  -  H  (2a  -  H)  (66) 

Similarly  I  the  equation  of  the  small  circle  of  the  planetary  cap 
having  height  H'  is: 

(x')2  +  (y')2  +  (z')2  -  a2 

z’  ■  a  -  H’  (67) 

x»2  ^  yi2  «  (2a  -  H’)  (o[ 

Using  Equation  (62)  we  can  write  Equations  (67)  and  (68)  in  terms  of 

the  (x,y,z)  coordinates# 

z’  •  X  cos  0c  +  z  sin  0c  “  a  -  H’  (69) 

x'2  +  y'2  _  0c  -  z  sin  0c)2  +  y  -  H'(2a  -  H')  (70) 


From  Equations  (65)  and  (69)  we  find: 

(a  -  H*)  -  (a  -  H)  cos  ^c 
^  Sin  0^ 

Using  this  result  and  Equation  (65),  Equation  (70)  gives: 


(71) 


2  2 
y  -  yp  “ 


H*(2a  -  H*)  sin^e^  -  [(a  -  H*)  cos  Qq  ~  (a  -  H)] 

8in^6(- 


(72) 


And,  of  course,  Equation  (65)  can  be  rewritten: 


z  ■  Zp  •  a  -  H 

Equation  (72)  requires  some  interpretation.  Figure  10  illustrates  a 
case  where  the  solar  and  planetary  caps  intersect#  It  is  possible 
that  the  caps  are  tangent  or  that  the  planetary  cap  lies  within  the 


(73) 
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solar  cap*  To  determine  whether  an  intersection  exists,  the  value 
of  yp^  of  Equation  (72)  may  be  used  as  a  discriminant* 


(1) 

If 

yp^ 

>  0, 

two  intersections  exist* 

(2) 

If 

yp^ 

-  0, 

the  caps  are  tangent* 

(3) 

If 

yp^ 

<  0, 

the  caps  do  not  intersect  and 
exposed  solar  angular  area  is 

Aex  "  0s  -  0p« 

the 

given  by 

For  the  great  circle  passing  through  the  intersection  of  the 
caps  z”  -  0,  z  -  zp,  X  ■  xp*  Hence,  Equation  (63)  gives: 

z”  -  ~  X  cos  6q  +  z  sin  Sq  "  0  (74) 

X  ^  /-rrx 

or  7  "  (75) 

thus  from  Equations  (71)  and  (73) 


(a  -  H*)  “  (a  -  H)  cos 
®G  “  (a  -  H)  Bln  Qq 


In  a  similar  manner  we  may  also  find  the  equation  for 
the  angle  between  the  center  of  the  planetary  cap  and  the  great 
circle^ 

tan  Oq’  -  “ 

everywhere  on  the  great  circle. 

Using  Equation  (62)  evaluated  at  x  ■  xp  and  z  ■  zp  we 
finally  obtain: 

(a  -  H')  cos  0Q  -  (a  -  H) 
tan  0Q  (a  -  H*)  sin  0q 


In  computing  Agx*  there  are  five  cases  of  Interest  and  these 
cases  are  Illustrated  In  Figure  12.  If  yp^  <  0  one  obtains  case 
V  where  Agx  ■  0^  -  0p4  If  yp^  >0,  however,  the  caps  Intersect; 
and  before  AgX  be  determined,  one  must  calculate  Ag  and 


(76) 


(77) 


(78) 
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Case  III.  Aex  -  08  “  0P  +  Ap  -  Ag  Case  IV.  AEX  ■  08  ”  Ag 
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ONE  HALF  OF  THE 
SPHERICAL  LUNE, 


X 


ONE  HALF  OF  THE 
SPHEPvICAL  LUNE,  Ap 

CIRCLE  OF  PROJECTED 
PLANETARY  CAP 


CIRCLE  OF 
PROJECTED 
SOLAR  CAP 


GREAT  CIRCLE 


Figure  13.  Area  of  Spherical  Lune, 


A7 


Ap,  where  As  is  the  angular  area  of  the  smaller  lune  formed  by  the 
solar  cap  and  the  great  circle,  and  Ap  is  the  angular  area  of  the 
smaller  lune  formed  by  the  planetary  cap  and  the  great  circle. 

Figure  13  illustrates  the  geometry  and  the  variables  used  to 
compute  Ag.  The  method  of  computing  Ap  is  identical  except  that 
all  computations  are  done  in  the  (x',  y’ ,  z’)  coordinate  system 
instead  of  the  (x,  y,  z)  coordinate  system. 

The  area  of  the  surface  ABC  of  Figure  13  is  given  by: 


area  of  ABC  « 


^0  ©bC'*’) 

ad4)  f  asin6cl0, 
0  6^(0 


To  get  the  angular  area,  Ag,  one  must  double  the  area  of  ABC  and 
divide  the  result  by  a^.  Thus, 


(79. 


2 


^0  W) 

sin  6  d0 


0 


(80) 


or 


where 


As 


[0a(4>)]  -  cos  [eb(<(>)] 


d(t> 


0 


(81) 


4)  is  the  angular  displacement  in  the  x-y  plane, 
positive  counterclockwise, 

6^(4))  is  the  angular  displacement  of  side  BC  (an 
arc  of  the  great  circle)  from  the  z-axis, 
positive  counterclockwise, 

6b(4>)  is  the  angular  displacement  of  side  AC  (an 
arc  of  the  small  circle  of  the  solar  cap) 
from  the  z-axis,  positive  counterclockwise. 

Before  Equation  (81)  can  be  evaluated,  expressions  for  cos  [6^(4>)] 
and  cos  must  be  found.  The  equation  of  the  small  circle  of 
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the  solar  cap  Is: 


x2  +  y2  .  H  (2a  -  H) . 


(82) 


Using  the  following  transformation: 

X  ■  a  sin  0  cos  (J) 
y  ■  a  sin  0  sin 
2  ■  a  cos  0 

Equation  (82)  becomes: 

a^  (sin^0  cos2(t>  +  sin^0  sin^$)  -  U  (2a  -  H) 


(R3) 


or 


(  1  -  cos^e)  -  H  (2a  -  H) 


(84) 


For  the  small  circle  of  the  solar  cap,  0  ■  Oj,,  and  Equation  (84) 
becomes 


(85) 


.2  -  a2 


(86) 


The  equation  of  the  great  circle  is 

(x")2  +  yi 

and  using  Equation  (53)  for  x"  we  obtain: 

(x  sin  Oc  +  2  cos  0g)^  +  «  a^.  (87) 

From  Equation  (75)  we  note  that  everywhere  on  the  great  circle 

z  tan  0G  *  x,  (88) 

Hence,  after  substitution  of  Equation  (88)  into  Equation  (87)  we  have 

(z  tan  Oq  sin  0G  +  z  cos  0g)^  +  y2  ■  (89) 

Inserting  the  polar  coordinate  relations  from  Equation  (83),  yields  after 
some  reduction: 

cos^0  (1  -  cos^0Q  8in^4j)  ■  cos^0g  cos^cfr  (90) 

or.  since  0  -  0a(4>)  along  the  great  circle 
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(91) 


cos  0r  COS 


COS  0a(4>) 


J ^  ~  cos^Sq  sin^(Ji 


Now  we  can  evaluate  Ag  of  Equation  (61)  using  the  value  of  cos 
and  cos  6^  above. 

fj 

0  r  cos  0p  cos  ((> 

. .  a  -  H 


As  ■  2 


r- 


cos^Gg  sin^^i 


d4> 


(92) 


Rearranging  we  have: 

diQ 


*5  -  2 


f  . ■  2  I  (1  -f)  d, 

J scc^Bq  -  sin^4)  o 


Letting  sin  ^  ■  x,  c^  ■  sec^^Q 


sin  (j>o 


Ag  ■  2 


dx 


.  2  ri  -  ill 


^0 


■  2  sin 


l"sin  4)0 

L 


sec  ©G 


-  2  (1  - 


-  2  sin“'  [sin  4*0  cos  Sg]  ~  2  <Po  f  1  “  ”) 


(93) 


From  Equation  (82)  and  Figure  13,  it  is  evident  that 


xp 


‘*’0  ■  /  2  2  "  (2a  -  H) 


/xp2  +  yp 


yp 


sin  ((>0  “ 


yp 


yxp2  +  yp2  (2a  -  H) 


,  yp  >  0 


(94) 


(95) 
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therefore , 


'J'O  ■  sin 


yp 

Al  (2a  -  H) 


0  1  ^0  1  f 


(96) 


Using  Equations  (93),  (95),  and  (96)  the  value  of  Ag  is  obtained. 
Notice  that  the  above  development  was  done  under  the  restriction, 

^  ^0  which  case  the  entire  lune  of  Intersection,  Ag,  is  in 

the  first  and  fourth  quadrants.  If  the  intersection  of  the  two  caps 
occurs  in  the  second  quadrant  (as  in  Case  I  of  Figure  12),  then  (po 
is  in  the  second  quadrant  and  the  entire  lune  of  intersection  lies  in 
the  second  and  third  quadrants.  Instead  of  recomputing  this  case, 

note  that  (see  Figure  13)  if  the  lune  of  intersection.  As,  lies  in 

the  second  and  third  quadrants  one  could  rotate  the  x-y  plane  about 
the  z-axis  by  180®,  resulting  in  the  problem  that  has  .lust  been  ana¬ 
lyzed.  Therefore,  Equations  (93),  (94),  and  (96)  give  the  correct 
value  of  the  angular  area  As  whether  <po  actually  lies  in  the  first 
or  second  quadrant.  Notice,  however,  that  if  ^  1,  4^0  .1^  then 
xp  *  0  and  cos  (pQ  ^  0;  but  if  ^  ^  then  xp  <  0  and 

cos  (pQ  ^  Thus,  cos  »Po  be  used  tc  discriminate  between  the 

two  cases  just  mentioned, 

A  completely  parallel  development  is  used  in  order  to  compute 
Ap.  The  major  difference  is  that  all  calculations  are  done  in  the 
(x’y’z’)  coordinate  system.  The  resulting  equations  are: 


Ap 


2  sin“^  [sin  4>o*  cos  0g*]  -  24)o*  [l 


(97) 


cos  (po* 


xp  cos  0^  —  (a  “  H)  sin  0^ 
Al'  (2a  -  II' ) 


(98) 


sin  4’ o' 


yp 

(2a  -  H') 


(90) 
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1 

■;)Q*  *  sin  ^ 

yp 

TT 

(2a  -  H' )  _ 

0  1  <Jo  1  - 

1 

wiiere  the  primed  quantities  are  related  to  the  planet  and  Of;*  is 
given  by  Equation  (78), 

Finally,  note  that  many  quantities  are  expressed  in  terms  of  , 
the  radius  of  the  great  sphere,  and  H  or  H*  the  depth  of  the  solar 
or  planetary  cap.  By  use  of  Equations  (60)  and  (61)  we  have 


U 

a 


27r 


H* 

a  ”  27r 

Thus  ail  quantities  may  be  expressed  in  terns  of  '.  <5  and  0p  which 

are  already  known.  This  is  most  easily  accomplished  by  settinp  n  ■=  1 

^  s 

(since  the  radius  of  the  sphere  is  arbitrary)  and  letting  11  »  ~ 

bp 

and  H’  -  “  . 

The  computation  of  the  exposed  angular  area  of  the  solar 

cap  is  shown  for  each  of  the  five  cases  of  Figure  12  in  Table  11.  Note 
that  the  discriminnts  used  to  determine  which  case  is  to  be  evaluated 
are  yp^t  cos  (;)o,  and  cos  4q** 

Finally  the  penumbral  illumination  factor  is  given  by  Equation 

(37): 

Aex 

P  'TT  • 


(lOO) 

(in;  'I 

(lOJ) 
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NUMLKICAL  INTECIIATION 


The  e(juacioiis  of  motion  for  a  space  vehicle  are  second-order 
dilfereniial  ecjUiitions  which  describe  the  accelerations  arising;  i  ror 
the  forces  acting  on  the  vehicle.  Accounting  for  the  primary  [jra\  i- 
cational  field  of  the  reference  body  and  four  types  of  perturbative 
acct lerat ions ,  the  equations  of  motion  (from  Equation  (2)), 


»  -  -  0  +  n  +  P2  +  P3  +  P.  (2) 

If  the  perturbations  are  considered  to  be  zero,  the  right  liand 

of  Equation  (2)  reduces  to  one  term  and  the  vehicle  will  follow’ 
a  Ecplerian  orbit  which  may  be  described  In  closed  form  in  terms  of 
true  or  eccentric  anomaly.  Usually,  however,  part  or  all  ol  tlir 
perturbations  are  included  and  Equation  (2)  must  be  numerically  inte- 
^jrated , 

There  arc  two  basic  methods  by  which  the  integration  of  Equation 
(2)  may  be  formulated  *  Encke's  method  and  Cowell’s  method.  If  Equa¬ 
tion  (2)  were  to  be  numerically  integrated  in  a  straight-forward 
manner,  the  integration  would  be  kno\>m  ns  Cowell's  method.  The  sim¬ 
plicity  of  this  method  is  offset  by  the  large  accelerations  which  must 
be  integrated.  As  a  consequence  of  the  acceleration  magnitudes,  small 
time  increments  have  to  be  used  in  the  Integration,  and  machine  round¬ 
off  error  accumulates  rapidly.  Independent  evaluations  at  many  com¬ 
panies  and  universities  have  shown  that  Cowell's  method  requires  more 

I.  u  U.11.  228-242] 

machine  time  than  other  perturbational  schemes.  See  Baker 

tor  a  further  discussion  of  Cowell's  advantages  and  disadvantages. 

Despite  its  drawbacks,  Cowell's  method  is  still  widely  used  and 

is  especially  suited  if  the  accelerations  due  to  perturbations  are 

as  large  or  larger  than  the  term  due  to  the  central  force  field  (e.r., 

during  reentry).  The  program  contains  both  methods  and  the  option  is 

left  to  the  user. 

Historically,  Encke's  method  is  older  than  Cowell's  although 
the  former  is  more  sophisticated.  Cowell's  method  requires  a  modern 
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high-speed  computer  to  be  practical  ^  whereas  r'n»  l:e’s  wa?;  drvclopc  tl 
for  hand  computation.  In  Lncke’s  method,  it  is  assumed  that  the  per¬ 
turbative  accelerations,  are  small  compared  to  the  reference 

body  acceleration.  Consequently,  when  drnr  accclern tions  arc  small, 
rhe  solution  of  Equation  (1)  is  a  good  approximation  to  the  true  orbit 
Under  these  conditions,  it  is  only  necessary  to  integrate  the  differ¬ 
ence  between  the  accelerations  on  the  two-hody  orbit  and  the  total 
accelerations  acting  on  the  vehicle.  The  equation^v  of  motion  then 
become  second-order  differential  equations  descrihinp  the  acceleration 
differences.  Let 

^  -  R  -  Rjjj 

where  R'pg  is  the  position  of  the  vehicle  in  ternvS  of  the  two-hody 
orbit.  Then, 


Equation  (lOA)  is  integrated  to  obtain  E,  and  These  quantities 

are  then  added  to  Rtb  and  RtB»  respectively,  to  obtain  the  instan- 
taneous  position  (R)  and  velocity  (R)  of  the  vehicle.  The  quan¬ 
tity  f  is  commonly  referred  to  as  the  ”Encke"  term. 


Encke’s  Method 

Rewriting  Equation  (104)  we  have. 


I 


Rtb  + 


The  above  equation  could  be  integrated  directly;  how^ever,  the 
Rr  R 

term  ^1  -  is  not  well  suited  for  the  nurerical  calculation 

3 

since  R^g/R^  is  very  close  to  1.  Hence  we  rewrite  this  term  bv 
setting 


(103) 


(104) 
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p  * 


-  TT 


-  i  - 


"TB 


R  •  R  “f'  (R  •  ^TB *  *  ^  *  ^TR  *  ^TB 

R  •  (Rtb  +  ^)  ~  Rfu  (RtB  +  c) 

- - - - -5 - 


3i  ,  -> 

Then,  (1  +  q)  ^  “  R^/R^B 


"TB 
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[(1+q)  ^  +  1]  1  -*-  Or  )'^^ 


^xpar.dinf;  c'.l  nurnerator 


R- 


^B 


-  1 


3q 


1  4- 


^  -3 

3q‘  ■*■  C 

- VT  • 

(1*^0  ^ 


In  bumnary,  LucKe*s  meihod  computes  the  ueviations  from 
orbit  by  integrating 


r3q  +  3q2  +  q3^,  ^ 

I - 57"  Rtb  - 

L 1  +  (]+c)  2  J 


V 

i’’* 


the  nominal 


where 


(2Ktb  + 

Rtb 


Determination  of  Keplerian  Orbit  Vectors 

Lncke’s  method  requires  the  calculation  of  the  position  and  velo- 

*  • 

city  vectors,  ^^TB»  respectively,  of  the  nominal  Keplerian 

orbit.  The  position  and  velocity  at  time  t  may  be  written  in  terms 
of  rne  position  and  velocity  at  an  earlier  time  Cq  follows: 

• 

RtbCc)  “  f  RTBCCq)  +  B  RTB(to) 

and 

•  • 

RTB(t)  “  ^  Rtb  (to)  8  Rtb  (to) 


wliure  f  and  g  are  explicit  functions  of  the  differential  eccentric 
anomaly  of  the  Keplerian  orbit.  The  program  uses  Herrick's  method  to 
determine  f,  g,  and  g,  and  then  computes  RTB(t)  and  RTB(t) 


(105) 

(106) 


(107) 

(108) 
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from  Equations  (107)  and  (108),  The  equations  for  computing  f,  f,  g, 
and  g  are  given  below;  their  derivations  may  be  found  in  Battin^^^ 
or  Baker^^°^. 


where 


where 


U  -  X'’  [ 


f  -  1  -  C/rg 

(109) 

Tq  "  |RTB(to)l 

(110) 

M  -  U 

8  " 

/y 

•  -  AT  s 

A  rg 

(111) 

g  -  1  -  C/a 

(112) 

/T”  (t-to)  -  ro  X  +  do  c  +  coU 

(113) 

r  .  JLi  +  Jtr  -  +  1 

'■2\  4!a  6!a^  8'a^  +  .  .  .  J 

(114) 

r  _  JLi  +  _  -Xi-  +  1 

^3!  5!a  7!a^  9!a^  •  •  • J 

(115) 

c  V  U 

S  "  X  -  ^ 

a 

(116) 

RtbC^o)  •  ®TB(to) 

“  J- 

/p 

(117) 

A  ■  ro  +  doS  +  coC 

(118) 

rg  v§ 

CO  -  -  1 

(119) 

i  -  -L  _  li 

(120) 

a  ro 

a  is  the  semi-major  axis  of  the  Keplerian  orbit, 
p  is  the  universal  gravitational  constant. 


and 
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Vo  •  I^TB(to)|. 


Solution  of  Equations  (109)  through  (112)  rsqulras  th«  dstsrul- 
nation  of  X.  Equation  (113)  la  solvad  for  X  by  writing 

F(X)  -  M  -  roX  -  dflC  -  cqU  •  0 

and 

r'(X)  -  -  ro  -  doC  -  CqU'  •  0. 

Salactlng  an  Initial  astinata  of  X  to  ba  X|,  tha  root  of 
F(X)  ■■  0  may  be  determined  by  a  Nawton-Raphaon  Itaratlva  procaaa, 

1  •  a  •  p 

•  Xi  -  *^^^iVy'(Xi) 

From  Equations  (114)  and  (115) 

C  -  X  -  U/a 
U'  •  C, 

Slnca  Equations  (114)  and  (115)  ara  aarlea  expanalona  in  pawara 
of  a  method  muat  be  found  to  Halt  tha  alra  of  thla  tava.  For 

an  elllpaa, 

X  -  (E  -  Eo)  /T 

where  £,  and  Eq  are  the  eccentric  aocaallaa  at  tlaea  tp  aad  to 

reapectlvely .  Hence p  by  updating  tha  epoch  to  at  fraquant  Intervale p 

the  difference  E  -  Eq  Is  kept  email.  By  slallar  raasonleg.  fre« 

X^ 

quent  updating  will  keep  small  for  hyperbolic  orbits. 

Cow£llJ^8__Metho^ 

As  described  beforsp  tha  general  equations  of  aetlon  ef  a  apace 
vehicle  are:  ^ 

^  I  ^1  (121) 

1-1 

In  Cowell's  method p  these  equatleas  are  Intagratad*  using  nuaerl-' 
cal  techniques*  to  obtain  the  Instaataaeeus  poaltlen  and  velocity  of 
tha  vehicle. 

St 


The  program  performs  the  integration  using  the  same  techniques  as  it 
does  for  the  Encke  method.  The  Runge-Kutta  starting  procedure  pro¬ 
vides  the  initial  data^  and  the  Nordseick  method  is  used  as  the  long¬ 
term  integration  procedure. 

Integration  Technique 

Equation  (2)  in  Cowell's  method  or  Equation  (105)  in  Encke 's 
method  must  be  numerically  integrated  by  the  program.  This  process 
is  divided  into  two  stages ^  a  starting  procedure  and  a  long-term  pro¬ 
cedure.  The  long-term  numerical  integration  procedure  requires  know¬ 
ledge  of  previous  data  points.  Thus,  the  starting  procedure  is  needed 
to  provide  the  initial  data  points  for  the  long-term  procedure. 

The  procedure  chosen  by  Sperry-Rand  for  long-term  integration  was 
based  on  their  experience  with  methods  used  in  the  ITEM  and  MINIVAR 
programs,  and  the  results  of  comparative  testing.  The  justification 
of  their  choice  is  presented  below  as  it  appears  in  the  original  docu¬ 
ment  on  the  SPACE  program^ 

**The  long-term  numerical  integration  procedure  presently  in  use  in 
the  ITEM  and  MINIVAR  programs  is  an  Adams  sixth-order  predictor  method 
(without  corrector)  for  second-order  differential  equations.  It  was 
desired,  however,  to  test  a  broader  class  of  procedures  before  deciding 
on  one  for  use  as  the  long-term  numerical  integration  procedure  to  be 
used  in  the  program.  Accordingly,  a  program  was  written  to  test  pre¬ 
dictor,  predictor-corrector,  and  iterated  predictor-corrector,  i.e., 
repeated  application  of  correctors,  techniques  of  various  orders  of 
approximation  and  with  or  without  modifiers. 

*^he  following  results  were  obtained: 

(a)  Modifiers  were  found  to  leave  the  error  unchanged.  There 
is,  therefore,  no  reason  to  use  them. 

(b)  As  time  interval  increases,  there  is  more  tendency  for  the 
solution  to  become  unstable.  Error  increases  with  a  large 
power  of  the  time  interval. 

(c)  As  the  degree  of  the  approximating  polynomial  increases,  a 
decrease  in  stability  is  noted.  Error  decreases  about  3:1 
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for  unit  increase  in  degree  of  the  approximating  polynomial, 

(d)  Predictor-only  methods  (degree  and  time  interval  held  fixed) 
are  about  30:1  less  accurate  than  predictor-corrector  meth¬ 
ods,  i.e.p  one  application  of  the  corrector  removes  97%  of 
the  error  in  the  predictor,  A  second  iteration  of  the  cor¬ 
rector  does  not  reduce  the  error,  but  does  Improve  stabil¬ 
ity  at  the  expense  of  a  2:1  increase  in  running  time.  The 
increase  in  running  time  is  intolerable;  therefore,  a  pre¬ 
dictor-corrector  method  will  be  used  with  no  iteration  of 
the  corrector, 

(e)  Either  a  fifth  or  sixth  degree  approximating  polynomial 
yields  a  good  compromise  between  accuracy  and  stability, 

**The  final  choice  of  a  long-term  integrating  procedure  involves 
considerations  other  than  only  accuracy  and  stability: 

(a)  The  ease  of  transforming  the  output  of  the  starting  proce¬ 
dure  to  the  form  of  input  starting  data  needed  by  the  long¬ 
term  procedure, 

(b)  Whether  the  long-term  procedure  can  easily  accomodate  a 
change  in  the  time  interval, 

(c)  Whether  the  long-term  procedure  can  easily  interpolate 
to  find  conditions  at  an  intermediate  time  at  which  data 
are  desired, 

**There  are  at  least  three  forms  in  which  the  Adams  long-uerm  pre¬ 
dictor-corrector  formulas  can  be  written.  The  only  difference  is  in 
mathematical  form;  therefore,  the  accuracy  and  stability  are  the  same 
for  all  three, 

*'The  first  form  is  the  conventional  one  in  terms  of  the  successive 
backward  differences;  the  second  is  in  terms  of  the  successive  values 
of  the  function.  The  third  is  due  to  Nordsieck  and  uses  the  succes¬ 
sive  higher  derivatives  of  the  approximating  polynomial.  Each  of  the 
three  forms  has  certain  advantages  and  certain  disadvantages  which 
will  be  discussed  now, 

**The  backward  difference  form  is  fairly  easy  to  start  but  inter¬ 
polation  is  somewhat  difficult,  and  it  is  virtually  impossible  to 
change  intervals  except  by  use  of  the  starting  procedure.  The  suc¬ 
cessive  value  form  of  the  method  is  trivial  to  start  up,  but  inter- 
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polation  involves  the  Lagrangian  interpolation  formulas,  and  changing 
intervals  is,  again,  almost  impossible*  The  inability  to  change  in¬ 
tervals  immediately  after  starting  causes,  as  in  the  present  ITEM  pro¬ 
gram,  a  situation  where  the  starting  solution  is  called  28  times,  but 
used  only  7  times* 

**The  Nordsieck  method  is  fairly  difficult  to  start,  but  very  amen¬ 
able  to  arbitrary  changes  of  time  intervals  and  to  interpolation  to 
intermediate  points.  Five  points  are  all  that  are  needed  to  start 
after  a  change  in  time  interval  (of  about  A;l).  Due  to  its  versatility, 
the  Nordsieck  method  of  degree  5  (called  m  ■  6  by  Nordsieck)  without 
iteration  and  without  choice  of  interval  is  used  in  the  program.” 

Starting  Procedure  (Runge-Kutta-Gill.  RKG.  Method) 

The  program  uses  the  Gill  modification  of  Runge-Kutta*- 
for  the  starting  procedure.  Runge-Kutta  methods  are  widely 
used  and  the  differences  between  Gill's  method  and  others  are  minor. 
Therefore,  only  the  formulas  as  given  by  Gill  and  not  the  deriva¬ 
tions  are  presented.  The  method  was  chosen  because  it  introduces  some 
simplicity  and  error  reduction  over  similar  methods. 

Given  the  differential  equation 


and  the  initial  values  xq  and  y(xo)i  y(xo+h)  is  calculated  by 


(122) 


where 


yo  ■  y(xo) 


(123) 
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The  above  procedure  could  be  used  directly  to  compute  y(xo  +  h) ^ 
y(xo  2h)  ^  j  etCt  However »  Gill  has  introduced  a  “bridging” 

technique  between  one  entry  and  the  next.  This  modification  increases 
the  accuracy  with  little  increase  in  complexity.  The  new  procedure^ 
may  be  siunmarized  as : 

Given  xq*  yo»  y*  (xi)  -  f(x^>  y(x^)) 

calculate  for  i  •  0^  1^  2^  3 

Ki  -  a^(y' (xi)  -  b^  qi) 

^1+1  -  X^  +  Ti  h 

yi+1  •  y(xi+i)  -  yi  +  h  Ki 


qi+i  ■  qi  +  3 

Ki  -  ci  y' (xi) 

1 

^0  ”  2 

ai  ■■  1  - 

bo  -  2 

/I 

CO  -7 

.-x./I 

1 

To  -  7 

Tj  ■■  0 

32  ■■  1  + 

/I  >=2-1 

C2  -  1  +  /I" 

1 

T2  -  7 

1 

33 

63-2 

1 

C3  -  2 

T3  -  0 

For  the 

first  entry  into  the 

procedure^  ^0  " 

and  the  proce- 

dure  is  the  same  as  given  in  Equations  (122)  and  (123).  For  subse¬ 
quent  entries  I  qg  is  set  equal  to  the  previous  q4 . 

It  should  be  noted  that  the  program  must  apply  the  method  to  the 
three  second-order  differential  equations  given  by  Equation  (105) 
(Encke’s  method) ^  or  by  Equation  (121)  (Cowell’s  method).  The  program 
treats  the  problem  as  one  of  six  first  order  differential  equations 
given  by: 

yi  “  fi(t»  yi»  . . ye)  i  ■■  i»  2 . 6 

where 


(12A) 
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yi (t)  -  x(t) 

y4(t)  -  x(t) 

y2(t)  ■  y(t) 

y5(t)  -  y(t) 

y3(t)  -  z(t) 

y6(t)  -  z(t) 

where  (x^  y,  z)  and  (x,  y,  z)  are  the  geocentric  position  and  velo¬ 
city  vectors  (in  Cowell’s  method)  or  the  perturbations  from  the  nomi¬ 
nal  orbit  (in  Encke’s  method). 

The  procedure  is  programmed  as  follows: 

1.  i*1^2^...^6  ^ 

2.  y  ■  (yi»  y2»  •••#  ye)  ■  (R(to)»  R(to)) 

3.  For  k  -  1,  2,  3,  4 

_•  4.  ^ 

A.  y  -  (yi.  ya.  ••••  ye)  (R»  R) 

B.  For  i  ■  1,  2,  , , , ,  6 

1.  K  ak  (yi  -  bk  •  qi) 

2.  yi  yi  +  K  •  h 

3.  qi  +  3  •  K  -  ck  •  yi 

c.  R  (yi,  ya.  ya) 

• 

R  (y4.  ys*  ye) 

t  t  +  Tk  •  h 

••  • 

call  DERIV,  i.e,,  cal.  ^  -  F(R,  R,  t) 

•  •• 

Hence  the  program  exits  with  R(to  +  h)  »  R(to  +  h) ,  and  R(to  +  h) . 
Also,  steps  1  and  2  are  only  executed  upon  the  first  entry  into  the 
procedure.  For  second  or  later  entries,  the  initial  value  of  is 

the  value  computed  during  the  last  entry  for  K  ■  A.  Clearly,  y 
already  contains  (R(t),  R(t))  and  need  not  be  reset# 

Transition  from  Starting  Procedure  to  Long-Term  Integration 

The  starting  procedure  yields  the  solutions  of  the  six  differen¬ 
tial  equations  and  their  rates  of  change  at  six  successive  times.  It 
is  necessary  to  transform  these  data  into  the  form  required  by  the 
Nordsieck  long-term  numerical  integration  procedure. 
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Bor  each  first  order  differential  equation,  the  Nordsieck  method 
requires  the  following  five  higher  derivatives  evaluated  at  t  ■  tgJ 


a(to) 

h  y(to) 
2! 

b(to) 

h2y(III) 

(to) 

3! 

c(to) 

h3  y(IV) 

(to) 

A! 

d(to) 

y(V) 

(to) 

5! 

e(to) 

h5y(VI) 

(to) 

6! 

J 

The  RKG  starting  method  provides  data  for  y(t)  and  y(t)  at 
the  six  time  intervals  up  to  and  including  to,  i«e«,  RKG  provides: 


y(t;o) 

y(to) 

y(to  “ 

h) 

ft 

o 

1 

h) 

y(to  “ 

2h) 

y(to  “ 

2h) 

y(to  - 

3h) 

y(to  “ 

3h) 

y(to  “ 

Ah) 

ft 

o 

1 

Ah) 

y(to  - 

5h) 

y(to  “ 

5h) 

The  required  values  for  a(to),  b(to)»  c(to)t  d(to)#  and  e(to) 
will  be  found  by  using  Lagrange's  Interpolation  formula  to  fit  a 
power  series  of  degree  five  to  the  y(t)  data  provided  by  the  RKG 
method.  The  power  series  will  then  be  successively  differentiated  to 
obtain  the  required  data.  Let 

X  -  t  -  to 

h 

and  let  primes  denote  derivatives  with  respect  to  x.  Therefore, 


(125) 
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f  (t)  -  h  y(t) 
y"(t)  -  y(t) 
y"'  (t)  -  h3 
y""(t)  -  h**  f 

y""’(t)  -  h5  y(VI)(t) 
y"""(t)  -  h® 


From  Lagrange’s  Interpolation  Formula, 


where 


Fo(x) 


Fl  (x) 


F2(x) 
FjCx) 
F4  (x) 


FsCx) 


5 

y(x)  -  I  Fi(x)  yi 
i“0 

(xfl)  (xf2)  (x+3)  (x+A)  (x^•5^ 

120 

-  x(x-h2)  (x+3)  (x+A)  (x+5) 

2A 

x(x+l) (x+3) (x+A) (x+5) 

12  \ 

-  x(x+l) (x+2) (x+A) (x+5) 

12 

x(x+l) (x+2) (x+3) (x+5) 

2A 

-  x(x+l) (x+2) (x+3) (x+A) 

120  -J 


and  y^  are  the  values  determined  by  the  RKG  method.  Multiplying 
out  the  factors  in  Equation  (128)  yields 


(126) 


(127) 


(128) 
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Fq  (x) 


Fl  (x) 


x^.  -f  13  x**  +  85  x\-f  225  +  274  x  +  120 

5  J.  lA  ^4 


1^ 

-  -  +  71  x3  4-  134  x"  +  120x) 

24 

S  4.  in7 


F2(x)  -  ^  ^  ^  ^  •^- 

^12  ■i'-  -f  >3  4-  78  x"  -f  Wli) 


F4(x) 


24 

120 


F5(x) 

From  Equations  (125)  and  (126) 

«(to)  "  ^  2'P'^ 

b(to)-^ 

c(to)-^ 


d(to) 


e(to) 


4! 

v'"'  (0) 

5! 


^Illl  f 


6t 


iOi 


Successively  differentiating  Equation  (127)  with  reepect  to  x, 
(using  Equation  (129)),  setting  x  ■  0,  end  substituting  into  Equa¬ 
tion  (130)  yields  the  following  in  matrix  notation; 


2a(to)~ 

'-24 

150 

-400 

600 

-600 

274~ 

y(to  - 

5h) 

3b(to) 

-50 

305 

-780 

1070 

-770 

225 

y(to  - 

4h) 

4c(to) 

.  -L.  . 

-35 

205 

-490 

590 

-355 

85 

y(to  - 

3h) 

5  d(to) 

120 

-10 

55 

-120 

130 

-  70 

15 

rt 

o 

1 

2h)| 

_6e(to)_ 

-  1 

5 

-  10 

10 

-  .5 

1 

rt 

o 

h) 

or 

o 

4J 

_ 1 

— 

(129) 


(130) 


66 


a(to) 

b(to) 

c(to) 

d(to) 

^e(to) 


1 

1440 


-144 

900 

-200 

1220 

-105 

615 

-  24 

132 

-  2 

10 

-2400 

3600 

-3120 

4280 

-1470 

1770 

-  288 

312 

-  20 

20 

-3600 

1644 

-3080 

900 

-1065 

255 

-  168 

36 

-  10 

2 

y(to 

-  5h) 

y(to 

-  4h) 

y(to 

-  3h) 

y(to 

-  2h) 

y(to 

-  h) 

y(to) 

(131) 


Equation  (131)  is  used  to  compute  the  sets  of  coefficients  [a(to)f 

T 

b(to)>  x(to)f  d(to)f  e(to)]  for  the  six  differential  equations  given 
by  Equation  (124). 


Long"-Term  Integration  (Nordsieck  Method) 

The  Nordsieck  method  is  used  to  continue  the  solution  of  Equation  (124) 


dyi 

-  fi(t,  . . ye)  i  “  1»  2, 

once  begun  by  the  RKG  method. 

Considering  the  single  equation 


dt 


f(t,  y) 


and  approximating  its  solution  by  a  polynomial  of  degree  five 

yP(to+h)  -  y(to)  +  h[y(to)  +  ^  y(to)  +  ^  y(^^^)(to) 


+  y(^'^)(to)  +  ^  y('^)(to)  +  ^  y('^^)(to)] 


5! 


6! 


(132) 


(133) 
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where 

h  is  the  integration  step  size, 
to  is  the  value  of  t  at  the  last  integration* 

Substituting  from  Equation  (125) 

y^(to+h)  -  y(to)  +  h[f(to,  y)  +  a(to)  +  b(to)  +  c(to)  +  d(to)  +  e(to)] 

Differentiating  Eqiaation  (133)  with  respect  to  t  =  t^  +  h  and  using 
Equation  (125) 

fP  -  f(to,  y)  +  2a(to)  +  3b(to)  +  4c(to)  +  5e(to) 


Having  obtained  a  predicted  value  of  f(to  +  b,  y(to  b))  a 
value  of  f(to  +  b)  ■  f(to  +  b,  y^)  is  computed  and  combined  with 
f^  to  give  tbe  Nordsieck  corrector, 


where 


Ki  h[f(to  +  h)  -  fP] 

^  19087  _ 

”  60A80  "  0.315591931 

The  values  of  a(to  +  h) ,  b(to  +  h) ,  c(to  +  h) ,  d(to  +  h)  and 


e(to  +  h)  are  then  computed  in  terms  of  their  values  at  t  ■  to 

by 


a(to  +  h)  -  a(to)  +  3b(to)  +  6c(to)  +  lOd(to)  +  15e(to) 

+  K2[f(to  +  h)  -  fP] 

b(to  +  h)  ■  b(to)  +  4c(to)  +  lOd(to)  +  20e  (to) 

+  K3[f(to  +  h)  -  fP] 

c(to  +  h)  -  c(to)  +  5d(to)  +  15e(to)  +  K4[f(to  +  h)  -  fP] 

d(to  +  h)  -  d(to)  +  6e(to)  +  K5[f(to  +  h)  -  fP] 

e(to  +  h)  -  e(to)  +  K6[f(to  +  h)  -  fP] 

where 

K2  -  -  1.141666667 


K3  -  I  -  0.625 

Ki.  -  ^  -  0,1770833333 


(134) 


(135) 


(136) 
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0.025 


Ks 


X. 

AO 

rr^  -  0.00138888889 


The  Nordsleck  method  may  now  be  summarized  as  follows.  For  each 

p  P 

entry,  y^  is  computed  by  Equation  (134),  f  by  Equation  (135),  and 

f(to  +  h)  -  f(to  +  h,  y^)  by  the  functional  relationship  implied  by 

Equation  (132).  Then  y(tQ  +  h)  is  found  by 

y(to  +  h)  •  yP  +  Kj  h(f(to  +  h)  -  f^) • 


Finally  the  coefficients  a(t)  through  e(t)  are  updated  by  Equa¬ 
tion  (136). 


The  integration  interval,  h,  is  readily  changed;  the  change 

being  accomplished  by  using  new  values  of  a(t),  b(t),  c(t),  d(t), 

and  e(t)«  These  new  values  are  obtained  from  the  following  equations: 

^n 

ho 


B 


an(t)  -  Bao(t) 
bn(t)  -  B2bo(t) 
Cn(t)  -  B^CgCt) 

dn(t)  -  B-do(t) 

en(t)  - 


where  the  subscripts  n  and  0  stand  for  new  and  old,  respectively. 
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COORDINATE  SYSTEMS  AND  TRANSFORMATIONS 

The  equations  of  motion  of  a  vehicle  are  given  by  Equation  (2)* 

It  is  necessary  to  express  the  vectors  used  in  this  equation  in  an 
inertial  coordinate  frame,  i*e,,  a  coordinate  system  in  which  inertial 
forces  due  to  the  system’s  acceleration  are  absent  or  at  least  negli¬ 
gible.  A  coordinate  system  rigidly  attached  to  the  Earth  is  inappro¬ 
priate,  since  such  a  system  experiences  noticeable  accelerations  due 
to  the  daily  rotation  of  the  Earth. 

The  inertial  frame  chosen  for  this  program  is  a  so-called  ’’mean 
equinox  of  base  date”  system  which  is  determined  by  the  vernal  equi¬ 
nox  of  0.0  hours  1  January  of  the  year  subsequent  to  the,  epoch  time 
(initial  input  time)  used  in  a  given  run.  This  particular  base  date 
system  has  been  chosen  as  a  basis  for  calculation  because  the  plane¬ 
tary,  lunar,  and  solar  coordinates  are  written  on  tapes  in  that  coor¬ 
dinate  system.  Rather  than  transform  the  tape  information,  the  vehi¬ 
cle  initial  conditions  and  the  oblateness  accelerations  are  transformed 
into  the  base  date  system. 

The  definition  of  the  base  date  system  and  various  auxiliary 
Earth  referenced  systems  is  given  below  along  with  the  transformations 
between  them. 

Terminology 

First,  consider  the  celestial  sphere  (Figure  14)  which  is  a 
sphere  of  infinite  radius  with  the  Earth  at  its  center  upon  which  the 
positions  of  stars,  planets,  the  Sun,  and  the  Moon  are  projected.  It 
is  sometimes  convenient  to  pick  the  center  of  the  sphere  to  be  at  an 
observer  or  the  center  of  the  Sun,  but  for  the  purpose  of  our  discus¬ 
sion  its  center  will  always  be  taken  to  be  at  the  geocenter.  The 
stars  may  be  taken  to  be  fixed  on  the  sphere  (their  small  proper 
motions  being  neglected)  while  the  Sun,  planets,  and  the  Moon  appear 
to  move  and  describe  certain  paths  along  its  inner  surface.  The  ap¬ 
parent  path  of  the  Sun  during  the  course  of  a  year  on  the  celestial 
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sphere  is  a  great  circle  called  the  ecliptic ,  The  north  celestial 
pole  is  the  point  where  the  extension  of  the  Earth’s  north  pole  inter¬ 
sects  the  celestial  sphere  and  the  south  celestial  pole  is  defined 
similarly.  The  great  circle  where  the  extension  of  the  Earth’s  equa¬ 
tor  intersects  the  celestial  sphere  is  called  the  celestial  equator. 

The  point  on  the  celestial  equator  where  the  apparent  path  of  the  sun 
crosses  into  the  northern  half  of  the  celestial  sphere  at  the  begin¬ 
ning  of  spring  is  called  the  vernal  equinox.  This  is  one  of  two 
points  where  the  ecliptic  and  equator  intersect.  In  the  following 
discussion  the  word  celestial  will  usually  be  dropped  and  we  will 
just  speak  of  the  north  pole  or  the  equator. 

The  equator  and  ecliptic  are  constantly  in  motion  due  to  the  ef¬ 
fects  of  nutation  and  precession.  It  should  be  pointed  out  that  this 
is  astronomical  precession  which  is  distinct  from  the  force  free  pre¬ 
cession  due  to  the  flattening  at  the  poles.  The  former  is  caused  by 
the  net  torque  on  the  equatorial  ’’bulges”  due  to  the  gravitational 
attraction  of  the  sun  and  moon.  The  torque  is  quite  small  so  that  the 
precession  is  extremely  slow  -  the  period  being  26,000  years  compared 
to  the  rotational  period  of  one  day.  The  total  applied  torque  is  not 
constant  in  time,  because  the  torques  of  the  Sun  and  Moon  have  slightly 
different  directions  to  the  ecliptic  and  vary  as  the  three  bodies  move 
around  each  other.  As  a  result,  there  are  irregularities  in  the  pre¬ 
cession,  commonly  designated  as  astronoinical  nutation.  The  astronomi¬ 
cal  nutation  is  not  to  be  confused  with  the  ’’true  nutation”  which  is 
present  even  in  the  force-free  precession  of  the  Earth’s  rotation 
axis.  In  the  subsequent  discussion  the  astronomical  precession  and 
nutation  will  simply  be  referred  to  as  precession  and  nutation. 

Precession  may  be  separated  into  luni-solar  precession,  and 
planetary  precession.  The  luni-solar  precession  causes  the  north 
pole  to  move  around  the  ecliptic  pole .  The  ecliptic  pole  is  the 
point  of  intersection  of  the  extension  of  the  normal  to  the  ecliptic 
plane  at  the  geocenter,  and  the  celestial  sphere.  Nutation  is  an 
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lA-746  D-B4 


NORTH  CELESTIAL  POLE 


ECLIPTIC 


CELESTIAL  EQUATOR 


Figure  14,  Celestial  Sphere 
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lA-748 


Figure  15.  Precession 
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Irregular  motion  of  the  pole  with  a  period  of  about  18.6  years.  These 
two  effects  combine  to  model  the  motion  of  the  true  equator.  The  ac^ 
tual  equator  at  any  time  Is  called  the  true  equator  of  date.  The  word 
date  here  refers  to  any  arbitrary  time.  The  mean  equator  of  data  la 
defined  by  the  location  of  the  equator  when  the  effects  of  nutation 
are  removed  (l.e.|  where  the  true  equator  of  date  would  be  If  there 
were  no  nutation)* 

The  gravitational  action  of  the  planets  on  the  Earth  causes  the 
plane  of  the  Earth's  orbit  to  slowly  precess.  This  motion  appears 
from  the  earth  as  a  slow  precession  of  the  ecliptic,  and  Is  called 
planetary  precession.  This  effect  causes  the  obliquity  (the  angle 
between  the  equator  and  the  ecliptic)  to  decrease  at  the  rate  of 
about  47**  a  century.  Therefore,  to  be  precise  we  refer  to  the  eclip¬ 
tic  of  a  specified  date. 

The  effect  of  precession  Is  Illustrated  In  Figure  15.  The  mean 
pole,  mean  equator,  and  ecliptic  at  time  to  labeled  base  date 

and  at  a  later  time,  t,  are  labeled  date.  Lunl-solar  precession 
causes  the  mean  equator  to  move  between  base  date  and  date,  while 
planetary  precession  moves  the  ecliptic.  Figure  16  Illustrates  how 
nutation  causes  the  true  equator  of  date  to  differ  from  the  mean 
equator  of  date. 

In  this  discussion  wa  will  use  several  coordinate  systems  which 
will  be  defined  In  terms  of  the  ^ifue  equator  and  ecliptic  and  mean 
equator  and  ecliptic  at  a  given  time.  The  true  vernal  equinox  (or 
true  equinox)  Is  defined  by  the  point  of  Intersection  of  the  true 
equator  and  ecliptic.  The  mean  vernal  equinox  (or  mean  equinox)  Is 
defined  by  the  point  of  Intersection  of  the  mean  equator  and  the 
ecliptic.  The  true  equinox  and  the  mean  equinox  are  points  on  the 
celestial  sphere  that  experience  a  slow  motion  and  hence  we  must 
specify  a  time  to  Indicate  their  exact  positions. 
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Figure  16.  Nutation 
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Coordinate  Systems 

The  five  coordinate  systems  used  by  the  program  for  generating 
trajectories  are: 

(1)  The  mean  equinox  of  base  date  system#  This  coordinate 
Vystem  fs  the  inert ial  system  emgTcyed  by  the  program 
and  employs  unit  basis  vectors  x,  y,  and  z  defined 
as  follows: 

X  is  a  unit  vector  directed  towards  the  mean  vernal 
equinox  of  base  date, 

z  is  a  unit  vector  normal  to  the  mean  equatorial 
plane  of  base  date,  positive  in  the  northern 
hemisphere. 

y  is  a  unit  vector  orthogonal  to  x  and  z, 

gqulnox  of  date  system.  This  system  employs 

unit  basis  vectors  xj,^,  y>^,  and  zm  defined  as 
follows: 

xj^j  is  a  unit  vector  directed  toward  the  mean 
vernal  equinox  of  date, 

zil  is  a  unit  vector  normal  to  the  mean  equatorial 
plane  of  date. 

yj.j  is  a  unit  vector  orthogonal  to  x^  and  zy, 

equinox  of  date  system.  This  system 

employs  unit  basis  vectors  xj,  yx»  and  zx 
defined  as  follows: 

XX  is  a  unit  vector  directed  toward  the  true 
vernal  equinox  of  date, 

z-p  is  a  unit  vector  normal  to  the  true  equatorial 
plane  of  date, 

y-p  is  a  unit  vector  orthogonal  to  x<p  and  zx« 

(4)  The  Greenwich  coordinate  system.  This  system 

employs  unit  basis  vectors  xGt  VGt  snd  zq 
defined  as  follows: 

XQ  is  a  unit  vector  directed  toward  the  inter¬ 
section  of  the  Greenwich  meridian  with  the 
true  equator  of  date, 

ZQ  is  a  unit  vector  normal  to  the  true  equator 
of  date. 
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yc  i«  a  unit  vector  orthogonal  to  xq  and  zq* 

(5)  The  topocentrlc  coordinate  system.  This  system  employs 

three  orthogonal  basis  vectors  Xg,  yg  and  Zg^ 
directed  toward  the  east>  north,  and  local  vertical, 
respectively,  at  a  point  on  the  Earth's  surface 
(considered  as  an  ellipsoid).  It  Is  only  used  for 
observation  calculations  and  its  details  are  dis¬ 
cussed  under  OBSERVATIONS* 

The  differences  between  the  above  coordinate  systems  are  due  to 
dynamical  effects  of  the  Earth's  motion  or  to  the  geometry  of  the 
Earth**  It  is  possible  to  express  a  3  ^  1  vector  of  position  or  velo¬ 
city  written  in  the  basis  of  one  coordinate  system  in  terms  of  another 

coordinate  basis  by  applying  3^3  orthogonal  matrix  transformations. 

— * 

Table  III  below  lists  the  transformations  used  by  the  program.  R  is  a 
3^1  vector  whose  subscript  Indicates  the  basis  in  which  the  vector 
is  expressed* 

Table  III. 

Transformation  Matrices* 


MATRIX 

EFFECT 

FROM 

TO 

EQUATION 

[G] 

Earth  Geometry 

ys 

XG 

yG 

ZG 

Rg  -  [G]  Rg 

[y] 

Right  Ascension  of 
Greenwich  from  True 

Equinox  of  Date 

XG 

yc 

ZG 

XT 

YT 

ZT 

[N] 

(Daily  Rotation) 

Nutation 

Xip 

H 

1 

Zt 

XM 

YM 

ZM 

-  [N]  Rt 

[P] 

Precession 

XM 

YM 

ZM 

X 

y 

z 

R  -  [P]  Rm 

The  derivations  of  [y]. 

[N] 

• 

and 

[P] 

are 

given  below. 

should  be  recognized  that  these  three  matrices  are  functions  of  time* 
Sidereal  Time 

When  computing  the  inertial  coordinates  of  a  point  or  station 
on  the  surface  of  the  rotating  Earth  at  a  particular  time,  a  relation¬ 
ship  must  be  found  between  the  Greenwich  meridian  (the  great  circle 
that  passes  through  the  true  north  pole  and  Greenwich)  and  the  vernal 
equinox  at  that  instant*  The  angle  measured  in  the  equatorial  plane 
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from  the  Greenwich  meridian  to  the  true  vernal  equinox  is  called  the 

apparent  Greenwich  sidereal  tiine« 

The  right  ascension  of  the  mean  equinox  measured  from  the  true 

equinox  is  called  the  equation  of  the  equinoxes  or  nutation  in  right 

ascension  and  is  usually  denoted  by  6a,  The  relationship  between 

the  apparent  sidereal  time  and  the  mean  sidereal  time  is: 

Mean  sidereal  time  -  Apparent  sidereal  time  -  6a, 

6a  is  always  measured  along  the  true  equator  from  the  true  equinox 

eastward.  The  apparent  and  mean  Greenwich  sidereal  times  are  tabu- 

[4]  V, 

lated  in  the  American  Ephemeris  and  Nautical  Almanac  for  0^  of 
each  day  of  the  year  along  with  the  equation  of  the  equinoxes,  6a. 
Each  entry  is  given  in  hours,  minutes  and  seconds  of  time  to  0®.001, 
where 

-  15°,  1™  -  15',  and  1®  -  15". 

Apparent  Greenwich  sidereal  time,  Oq,  since  it  is  measured 
from  the  true  equinox,  is  affected  by  nutation  and  hence  changes  at 
an  irregular  rate  (see  Figure  17).  Since  mean  sidereal  time, 
is  measured  from  the  mean  equinox,  which  is  affected  only  by  preces¬ 
sion,  it  increases  at  almost  a  constant  rate  with  only  a  small  accel¬ 
eration  and  may  be  calculated  at  0^  of  any  day  by  the  formula 
0G  -  6^^  38®  45!836  +  8 ,640,184!542T  +  0?0929t2 

where  T  is  the  number  of  Julian  centuries  (36,525  days)  from  noon 
0  January,  1900  (Julian  day  2,415,020.). 

An  equivalent  formula  in  degrees  is 

eG(d)  -  100;0755426  +  0:985647346  d  +  2:9015  x  10"^ ^  d^ 

where  d  is  the  integer  number  of  days  past  O^'  1  January,  1950, 

Equations  (137)  and  (138)  may  be  reformulated  as  a  function  of  the 
time  past  any  conveniently  chosen  epoch  resulting  in  an  equation  of 
the  same  form  with  different  coefficient.. 

Evaluating  the  above  formula  for  6G(d  +  1)  and  6G(d),  we 


(137) 


(138) 
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can  find  the  angle  through  which  the  Greenwich  meridian  rotates  in  a 
day  with  respect  to  the  moving  mean  equinox# 

0G(d  +  1)  -  eG(d)  +  360?  -  360!985647346  +  5?  x  d. 

in  which  the  additive  term  2?9015  x  10”^^  has  been  neglected*  Di- 
viding  Eqixation  (139)  by  the  number  of  seconds  in  a  day  (86,A00)  and 
neglecting  the  final  terra  whose  contribution  is  insignificant  over  a 
period  of  only  a  few  years  (  <  25  yrs#),  we  have  cjg,  the  rotational 
rate  of  the  Earth  with  respect  to  the  moving  mean  equinox. 

Us  -  0!0041780  74622  sec."l 

Therefore,  to  find  the  mean  Greenwich  sidereal  time  at  time 
d  +  T  we  use 

0G(d  +  t)  -  0G(d)  +  Ws  T 

where  d  represents  the  integer  number  of  days  past  0^  1  January, 

1950,  and  t  the  number  of  seconds  that  have  elapsed  since  0^  of 
the  day. 

The  apparent  Greenwich  sidereal  time  6q  can  be  found  at  time 

t  by 

0c(t)  -  0G(t)  +  6a. 

Now,  (Ug  will  be  referred  to  as  the  sidereal  rotational  rate 
of  the  Earth.  The  word  sidereal  is  used  since  the  rate  is  with  res¬ 
pect  to  the  mean  equinox,  which  is  the  reference  point  used  for 
sidereal  time. 

Transformations 

The  effects  of  nutation  and  precession  are  fairly  small  and 
over  a  period  of  a  few  days  the  errors  Introduced  by  neglecting  these 
effects  Is  on  the  order  of  ten  meters  at  most.  Therefore,  for  short 
trajectory  determination  runs  the  program  has  the  option  of  not  In¬ 
cluding  them  by  setting  the  [P]  and  [N]  matrices  equal  to  the 
Identity  matrix,  [I],  The  [y]  matrix  which  includes  the  effects 
of  dally  rotation  must  be  Included,  however. 


(139) 


(140) 


(141) 


(142) 
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Figure  17.  Sidereal  Time 
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.  Figure  18.  Formation  of  Nutation  Matrix 
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Referring  to  Figure  17  the  transformation  from  the  Greenwich  co¬ 
ordinate  system  to  the  true  equinox  of  date  system  consists  of  a  ro¬ 
tation  along  the  true  equator  through  the  angle  6G(t)  where  0G(t) 
is  given  by  Equation  (142).  Thus  the  [y]  matrix  is  given  by: 


[y] 


cos  0G(t)  -  sin  0G(t)  0 

sin  0G(t)  cos  0G(t)  0 

0  0  1 


(1A3) 


If  the  effects  of  nutation  and  precession  are  to  be  neglected  then 
replaces  ©g^^^  in  Equation  (143),  and  the  reference  system 
of  the  computation  may  be  considered  as  the  system  defirted  by  the 
mean  equinox  of  date  (i«e«,  the  effects  of  nutation  and  precession 
are  neglected) , 

Next,  an  expression  for  the  nutation  matrix  [N]  is  found. 

In  Figure  18,  y  is  the  true  equinox  of  date  and  y  is  the 
mean  equinox  of  date.  The  right  ascension  of  the  mean  vernal  equinox 
referred  to  the  true  equinox  is  6a,  the  equation  of  the  equinoxes. 
Notice  that  6a  is  less  than  zero  as  it  is  drawn  in  Figure  18,  since 
right  ascension  is  measured  positively  toward  the  east.  The  nutation 
in  longitude,  ^  is  the  longitude  of  the  mean  equinox  measured 
from  the  true  equinox.  Celestial  longitudes  are  measured  along  the 
ecliptic,  and  since  the  positive  direction  is  east,  is  also  neg¬ 

ative  in  Figure  18.  The  mean  obliquity,  e,  is  the  inclination  of 
the  ecliptic  to  the  mean  equator,  and  e,  the  true  obliquity,  is  the 
inclination  of  the  ecliptic  to  the  true  equator.  The  nutation  in 
obliquity  is  6e,  where 

6e  ■  e  -  e. 


From  inspection  of  Figure  18,  and  by  using  Napier's  rules  for 
right  spherical  triangles, 

sin(90®  -  e)  ■  tan  6a  tan(90®  -  AiJ^) 
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or 


tan  6a  -  tan  Aij/  cos  e 


Using  the  series  expansion  for  the  tangent  function  we  have 


6a  + 

Now, 

fore 


(6a)3  2  C6a)5, 

3  15  ^  * 

I  6a  I ,  and  |  6iJj  | 

,  Che  approximation 


,  .  .  eos  £  (1*  +  iMll  +  +  .  .  . 

are  less  than  li  x  lO""**  radians,  There- 


6a  -  A\)/  cos  c 


should  never  produce  an  error  greater  than  1  ^  10"^^  radians  or 
2  X  10  ^  seconds  of  arc* 

The  nutation  matrix  [N],  is  now  derived  which  transforms  co¬ 
ordinates  referred  to  the  true  equinox  of  date  into  coordinates  re¬ 
ferred  to  the  mean  equinox  of  date*  This  transformation  consists  of 
three  rotation  matrices*  Using  Figure  18,  we  note  that  the  first, 
[A],  is  a  rotation  about  the  x-axis  through  the  angle  e*  The  sec¬ 
ond,  [B],  is  a  rotation  about  the  z*-axis  through  AiJ;,  and  the 
third,  [C],  is  a  rotation  about  the  x”-axis  through  -  e. 
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0 


0 

•  In  c 
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[N]  -  [C][B][A] 


[C]  - 


10  0 
0  cos  G  -  sin  e 

0  sin  G  cos  G 


cos 

-  cos  G  sin 

-  sin  G  sin 


sin  cos  G 

cos  G  cos  cos  G 
+  sin  G  sin  g 

sin  G  cos  A^/  cos  g 
-  cos  G  sin  G 


sin  Ail^  sin  g 

cos  G  cos  A\|^  sin  g 
-  sin  G  cos  G 

sin  G  cos  A\|^  sin  g 
cos  G  cos  G 


(14A) 


Since  |6\p|  <  10”^  and  |g  -  g|  <  10“^  the  nutation  matrix  is 
often  approximated  by  the  following  matrix: 


[N'] 


1 

-  A\^  cos  G 

-  A\^  sin  G 


A\|/  cos  G 
1 

-  6g 


A\|;  sin  G 

6g 

1 


The  error  in  the  elements  in  [N*]  should  be  less  than  one  unit 
in  the  eighth  significant  figure^^^* 

Ref erence  [ l4]  gives  formulas  for  A\^  and  6g  which  depend  on 
the  mean  longitude  of  the  Sun,  the  mean  longitude  of  the  perigee  of 
the  Sun,  the  mean  longitude  of  the  Moon,  the  mean  longitude  of  the 
ascending  node  of  the  Moonp  and  the  mean  longitude  of  the  lunar  peri¬ 
gee.  There  are  69  terms  in  the  expression  for  A\^  and  AO  terms  in 
6g.  Truncated  expressions  for  and  6e  are  given  in  Reference  [l]. 

Finally,  an  expression  for  [P]>  the  precession  matrix  is  de¬ 
rived  below.  Figure  19  illustrates  the  effect  of  precession  over  the 
period  of  time  t  -  to*  In  this  figure,  90®  -  to  right  ascen¬ 

sion  of  the  ascending  node  of  the  mean  equator  at  time  t  on  the 
mean  equator  at  to  measured  from  the  mean  equinox  of  to;  90®  -f  z 
is  the  right  ascension  of  the  node  measured  from  the  mean  equinox  of 
t;  and  0  is  the  inclination  of  the  mean  equator  at  time  t  to  the 
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Figure  19.  Formation  of  Precession  Matrix 
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mean  equator  of  to* 

The  precession  matrix^  [?]»  is  now  derived  which  transforms 
coordinates  referred  to  Y(t)  to  coordinates  referred  to  Y(to)* 

[P]  is  a  result  of  three  rotation  matrices.  From  Figure  19  we 
see  that  the  firsts  [A’]^  is  a  rotation  of  (90®  +  z)  about  the  z- 
axis.  The  second ^  is  s  rotation  about  x  through  -0  and 


the 

third,  [C'l. 

is  n. 

rotation  about 

z 

through  (Co 

-  90“). 

[P]  -  [A’] 

[BM 

[C], 

where 

COS 

VO 

O 

0 

+  z) 

sin 

(90*  +  z) 

0 

-  sin  z 

cos  z 

[A’] 

- 

-  sin 

(90“ 

+  z) 

cos 

(90“  +  z) 

0 

- 

-  cos  z 

-  sin  z 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

[B’l  - 

0 

cos 

(-0) 

sin  (-0) 

- 

0 

cos  0 

-  sin  0 

0 

-  sin 

(-0) 

cos  (-0) 

0 

sin  0 

cos  0 

[C] 


^  cos  (i;o  -  90®) 
-  sin  (;o  -  90®) 
0 


sin  (;o  -  90®)  0 

cos  (^0  "  90®)  0 

0  1 


sin  ^0  "  CO®  ^0  0 

cos  ^0  ^0  0 

0  0  1 


Hence, 


sin  z 

sin  ;o 

COS  z 

sin 

40 

sin 

0 

cos 

40 

+  cos 

z  cos  0 

cos 

40 

+  sin 

z  cos  0  cos  ^0 

- 

-  sin  z 

cos  ^0 

cos  z 

cos 

40 

-  sin 

0 

sin 

40 

-  cos 

z  cos  0 

sin 

40 

-  sin  z 

cos 

0  sin  Co 

-  cos  z 

sin  0 

-  sin  z 

sin 

0 

cos 

0 

Ref erence  [  14]  gives  the  following  numerical  expression  for  ^o» 
and  0  for  the  case  when  coordinates  are  precessed  to  a  later  epoch. 
Initial  epoch,  to  ■  1900,0  +  Tq 
Final  epoch,  t  ■  1900,0  +  tq  +  t 

Co  •  (2304. "250  +  l."396  Tq)  t  +  0."302  +  0."018 

z  -  Co  +  0."791 

0  -  (2004. "682  -  0."853  To)  t  -  0."426  -  0."042 


(1A5) 
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where  tq  and  t  are  in  units  of  tropical  centuries. 

When  the  initial  epoch  is  1900.0  +  Tq  +  t  and  the  final  epoch  is 
1900.0  +  Tq  (i.e.,  we  are  precessing  a  vector  to  an  earlier  epoch) 
one  vhould  replace 
Co  by  -  2 
z  by  -  Co 
0  by  -  0 

in  the  above  matrices  [A*],  [B*],  [C*],  and  [P]. 

Combining  the  results  of  Equations  (143),  (144),  and  (145)  with 
the  relations  in  Table  III,  all  necessary  coordinate  transformations  may 
be  performed  by  the  program. 

OBSERVATIONS 

The  SPACE-A  program  has  a  provision  for  a  total  of  25  observation 
types  of  which  a  number  of  ground-based  observations  are  specified  and 
described  here.  The  specified  observations  include: 


azimuth,  A 

elevation,  E 

round-trip  range,  2p 

. 

one-way  range  rate,  p 

hour  angle,  H 

declination,  D 

i-direction  cosine,  i 

m-direction  cosine,  m 

X-angle,  X 

Y-angle,  Y 

range  equivalent  time.  At 

range-rate  equivalent  time.  At* 


The  position  and  velocity  of  a  vehicle  is  available  at  all  times 
from  the  trajectory  generation  portion  of  the  program.  The  station 
position  and  velocity  is  cpmputed  from  input  data  concerning  its  geo¬ 
detic  coordinates.  Therefore,  one  may  consider  the  relative  geometry 
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of  the  vehicle  and  station  to  obtain  the  observations. 


The  program  has  the  option  to  Include  the  effects  of  refraction 
upon  ground-based  observations.  A  simple  model  of  the  Earth* s  atmos¬ 
pheric  effects  upon  observations  Is  employed  In  which  the  Index  of  re¬ 
fraction  varies  only  with  altitude.  Neglecting  the  effects  of  refrac¬ 
tion  causes  the  predicted  elevation  angle  to  be  less  than  the  observed 
elevation.  Errors  are  also  made  In  the  computation  of  range  and  range 
rate  related  observations. 

The  method  used  to  Include  the  effects  of  refraction  Is  to  com¬ 
pute  range  and  range  rate  from  the  relative  geometry  of  the  station 
and  the  actual  vehicle  position  and  velocity  and  then  to.  apply  correc¬ 
tions  based  on  the  refraction  model  In  order  to  obtain  the  range  and 
range  rate  related  observations.  In  the  case  of  angle  related  mea¬ 
surements  a  slightly  different  approach  Is  taken.  A  vector  pointing 
from  the  station  to  the  computed  vehicle  position  Is  constructed. 

Then  the  refraction  model  Is  used  to  obtain  a  new  vector  pointing  from 
the  station  to  the  observed  vehicle  position.  All  angle  observations 
are  then  computed  directly  or  Indirectly  using  this  new  vector. 

Coordinates  Used  In  Observation  Calculations 

The  Input  coordinates  of  the  station  (the  station  Is  assumed  to 
be  on  an  ellipsoidal  Earth)  are  Initially  given  as  geodetic  longitude, 
(Xq),  geodetic  latitude  $  and  height.  Given  this  Information 

the  position  and  velocity  of  the  station  In  the  true  system  of  date 
coordinate  system  Is  given  by: 

Xj  ■  [c  +  h]  cos  cos  (Xq  +  y) 

yx  -  [c  +  h]  cos  <|)G  sin  (Xq  +  y) 

2x  ■  [c(l-c^)  +  h]  sin  4>q  ^ 

XX  -  -  Wg  yx 

yx  -  uje  Xx 

2X  ■  0 


(146) 
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c  ■ 


:=a= 

sin^  (^G 


(147) 


e2  -  2f  -  f2 

where 


(148) 


yT»  Z7  are  the  true  system  of  date  coordinates  of 
the  station, 

IGi  ^  geodetic  coordinates  of  tha  station, 

Y  Is  the  right  ascension  of  the  Greenwich 
meridian  In  the  true  system  of  date, 

Re  Is  the  equatorial  radius  of  the  Earth, 

f  Is  the  flattening  constant  of  the  Earth- 


- -  - - - - -  ‘-298. 30^ 

Using  the  precession  matrix,  [P]  and  the  nutation  matrix,  [N],  we 

can  express  the  station  coordinates  in  the  inertial  (base  date)  coor¬ 
dinate  system.  From  the  position  and  velocity  vectors  of  the  station 
and  the  vehicle  one  can  obtain  the  vectors  of  position  and  velocity  of 
the  vehicle  with  respect  to  the  station  expressed  in  the  inertial  sys¬ 
tem.  Thus 

P  -  Pg  -  R 

UL-i 


(149) 

(150) 


where  p,  p  are  vectors  of  the  vehicle  with  respect  to  the 
station, 

Psf  Pg  are  the  vectors  of  the  station,  and 
R,  ^  are  the  vectors  of  the  vehicle. 

Note  that  all  of  these  vectors  are  expressed  in  the  inertial  coordi¬ 
nate  system. 

For  observation  calculations  it  is  convenient  to  define  a  station 
oriented  coordinate  system  which  has  a  basis  consisting  of  three  unit 
orthogonal  vectors  pointed  toward  the  east,  north  and  local  vertical. 
Using  station  oriented  coordinates  the  unit  vectors  are: 
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Xg  -  [10 

0]'^, 

east  vector 

ys  "  [0  1 

0]'^, 

north  vector 

-  [0  0 

up  vector 

Expressed  in  terms  of  the  inertial  base  date  system  these  vectors  arc; 


e  -  [T]  Xg 

east  vector 

(151) 

n  -  [T]  yg 

north  vector 

(152) 

u  -  [T]  zs 

up  vector 

(153) 

and 

[T]  -  [P][N][Y][G][a] 

(154) 

where  all  vectors  are  3x1  vectors. 


[P].[N],  and 


[a]  is  a  3  X  3  rotation  matrix  to  account  for 
misalignments  between  the  station  oriented 
coordinates  and  the  true  topocentric  east^ 
north,  and  vertical  coordinates, 

[G]  is  a  3  X  3  matrix  transforming  coordinates 
expressed  in  the  topocentric  system  into 
coordinates  of  the  Greenwich  system, 

[y]  are  the  3x3  matrices  of  precession,  nuta¬ 
tion,  and  right  ascension  of  the  Greenwich 
meridian  which  transform  the  coordinates  into 
the  inertial  base  date  system. 


[G] 


sin  Xq  -  sin  (^q  cos  Xq 

cos  Xq  -  sin  4>G  sin  Xq 

0  cos  <i>Q 


cos  4>g  ^G 
cos  (J)G  sin  Xq 
sin  (J)G 


(155) 


where  Xq  and  <(>g  are,  respectively,  geodetic  longitude  and  geo¬ 
detic  latitude  of  the  station. 

The  final  result  of  the  above  transformations  is  that  the  station 
oriented  basis  vectors  (e,  n,  and  u)  are  expressed  in  the  inertial 
base  date  coordinate  system. 


Computation  of  Observation  Types 

In  the  following  discussion  all  computations  are  made  using  vec¬ 
tors  expressed  in  the  inertial  (base  date)  coordinate  system. 

Azimuth  and  elevation  measurements  are  illustrated  in  Figure  20. 
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Azimuth,  A,  is  measured  easterly  from  station  north  from  0  ^  A  ^  2it. 
Elevation,  E,  Is  measured  from  the  stations  horizontal  plane  upward 

TT  TT  -► 

with  a  range  ”  “  i  E  i  “•  If  P  1*  the  station  to  vehicle  vector 
the  azimuth  and  elevation  are  given  by: 


A  - 


E  ■  tan 


-1 


/ 


0  *  U 


(p  •  e)2  +  (p  •  n)- 


If  refraction  Is  to  be  Included  In  the  computations,  the  vector,  p. 
Is  replaced  by  the  vector  p' ,  where  p'  Is  the  vector  from  the 
station  to  the  observed  vehicle  position  as  described  below. 


(156) 


(157) 


Figure  20.  Azimuth  and  Elevation 


The  round-trip  range,  2p,  Is  twice  the  distance  from  the  station 
to  the  vehicle  and  Its  value  Is  given  by: 

2p  -  2|p|  +  2Ap  (158) 

If  the  effects  of  refraction  are  not  Included  In  the  calculation 
2Ap  ■  0.  If  the  refraction  effects  are  Included  Lp  ia  a  correction 


92 


term  the  value  of  which  is  given  by  Equation  (202). 

t 

The  one">vay  range-^^ate^  pp  is  equal  to  the  magnitude  of  the 
time  derivative  of  the  vector  from  the  station  to  the  vehicle  and  its 
value  is  given  by;  • 


P 


4*  Ap 


(159) 


If  the  effects  of  refraction  are  not  included  Ap  ■  0.  If  the  effects 
of  refraction  are  included  Ap  is  a  correction  term  the  value  of 
which  is  given  by  Equation  (206) • 

The  hour  angle,  H,  is  the  angle  between  the  station  meridian 
and  the  projection  on  the  true  equator  of  the  station  to  vehicle  vec¬ 
tor  measured  in  the  true  equatorial  plane.  It  is  measured  westward 
from  -  TT  ^  H  ^  TT,  The  declination,  D,  is  the  angle  made  with  the 
true  equatorial  plane  by  the  station  to  vehicle  vector.  Declination 
is  measured  positive  in  the  northern  hemisphere  and  has  a  range 
-  Y  ^  D  ^  Hour  angle  and  declination  are  illustrated  in  Figure  21 
which  shows  the  station  at  the  center  of  a  celestial  sphere  and  the 
projected  true  equator  on  the  celestial  sphere.  The  hour  angle  and 
declination  may  be  derived  using  spherical  trigonometric  relations  in 
terms  of  geodetic  latitude,  azimuth  and  elevation  (4),  A,  E) ,  The 
values  are  given  by: 


H 


tan*^ 


sin  A 

cos  A  sin  -  cos  4)^  tan  E 


(160) 


D  -  sin"^  [sin  4>G  ^  cos  E  cos  A] 


(161) 


The  values  of  A  and  £  used  above  are  given  by  Equations  (156)  and 
(157). 

The  jZ>-direction  cosine,  £,  and  m-direction  cosine,  m,  are 
shown  in  Figure  22,  The  ^,-direction  cosine  is  the  cosine  of  the  angle 
between  tne  station  to  vehicle  vector  p,  and  the  station  east  vector, 
e.  The  m-direction  cosine  is  the  cosine  of  the  angle  between  the  sta- 
tion  to  vehicle  vector  p,  and  the  station's  north  vector,  n.  The 
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values  of  I  and  m  are  given  by: 


(162) 


m  ■ 


(163) 


If  refraction  effects  are  included  p’  replaces  p  in  the  above 
equations# 

^"^Pglg  and  Y-angle  measurements  are  illustrated  in  Figure  23. 
The  Y-aingle  is  the  angle  between  the  station  to  vehicle  vector^  p^ 
and  the  perpendicular  projection  of  this  vector  on  the  station’s  east- 
vertical  plane.  It  is  positively  measured  toward  the  north,  negatively 
toward  the  south  with  limits  "  ^  ^  ^  X-angle  is  measured 

between  the  vertical  vector,  u,  and  the  perpendicular  projection  of 
the  station  to  vehicle  vector  onto  the  station’s  east-vertical  plane. 

It  is  measured  from  the  positive  vertical  eastward  with  limit 
-  X  ^  TT,  The  value  of  X  and  Y  are  given  by: 


X  - 


Y  -  tan 


-1 


/ 


P  »  n 


(p  •  e)2  +  (p  •  u)' 


If  refraction  effects  are  included  p’  replaces  p  in  the  above 
equations. 

The  range  time  equivalent.  At,  and  the  range-rate  time  equiv¬ 
alents  At’,  are  included  since  the  raw  data  from  typical  tracking 
systems  are  the  time  between  a  transmitted  and  received  signal  for 
range,  and  the  time  to  count  a  given  number  of  Doppler  cycles  for 
range-rate.  In  most  systems,  these  quantities  are  first  converted  to 
range  and  range-rate  units.  However,  it  may  be  found  useful  in  some 
cases  to  use  the  raw  data  equivalents.  The  values  of  At  and  At’ 
are  given  by: 

At  ■ 

c 


(16A) 


(165) 


9A 


(166) 


Zenith 


Figure  21.  Geometry  Illustrating  Hour  Angle  ard  Declination 


Figure  22.  -cosine  and  m-cosine  Figure  23.  X-Angle  and  Y-Angle 
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At'  - 


(167) 


c 


where 

c  is  the  velocity  of  light, 

[pil  and  IP2I  are  the  ranges  at  the  beginning  and  end 
of  the  measurements ,  respectively • 

Specific  provision  is  made  in  the  program  for  computing  dopplcr 
shift  frequency  measurements.  The  equations  for  these  observations 
are  not  included,  since  the  method  of  computation  can  change  depending 
on  the  specific  hardware  in  use  for  a  particular  mission. 

Effects  of  Refraction 

Because  of  variations  in  the  refractive  index  of  the  atmosphere 
the  propagation  path  of  a  tracking  beam  is  subject  to  refractive 
bending. 

Two  models  for  the  variations  of  refractive  index  are  included 
in  the  program:  one  is  a  model  for  the  troposphere,  the  other  is  a 
model  for  ionosphere.  A  numerical  method  due  to  Weisbrod  and 
Anderson^ is  employed  in  which  the  effects  of  refractive  bending 
are  determined  by  numerically  integrating  over  the  total  propagation 
path,  the  index  of  refraction  at  each  point  being  determined  by  the 
assumed  model.  The  vector  from  the  station  to  the  current  calculated 
position  of  the  vehicle,  p,  is  available  from  the  program.  Using 
numerical  methods  a  correction  to  the  elevation  angle  is  found  and  a 
new  vector,  p',  from  the  station  toward  the  observed  position  of 
the  vehicle  is  determined. 

In  addition  to  angle  corrections  due  to  refractive  bending,  the 
effects  on  range  measurements  due  to  signal  retardation,  and  the  ef¬ 
fects  on  range— rate  measurements  due  to  refractive  bending  are  also 
found • 

Index  of  Refraction  Models 

In  order  to  simplify  computational  problems,  atmospheric  models 
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representative  of  average  conditions  are  employed  in  which  the  fol¬ 
lowing  assumptions  are  made) 

(1)  The  gradient  of  the  index  of  refraction  varies  only 
with  altitude. 

(2)  The  index  of  refraction  profile  is  approximated  by  a 
number  of  linear  segments  in  which  the  length  of  each 
segment  is  very  small  compared  to  the  Earth's  radius, 

(3)  The  troposphere  extends  to  40  kilometers, 

(4)  The  region  between  the  end  of  the  troposphere  and  the 
beginning  of  the  ionosphere  is  assumed  to  have  zero 
refractivity  (and  therefore  no  bending  or  signal  re¬ 
tardation  occurs), 

(5)  The  ionosphere  lies  between  a  height  ho  (input  data) 
and  2^000  kilometers, 

(6)  The  index  of  refraction  is  zero  beyond  2,000  kilometers. 

In  the  tropospheric  model,  refractivity  (N)  is  assumed  to  decay 

exponentially,  with  the  ground  index  of  refraction  and  the  scale 
height  as  parameters.  The  equation  for  the  tropospheric  model  is  as 
follows : 

N  -  No  -  (n  -  1)10^ 

where 

Nq  is  the  refractivity  at  sea  level,  an  input  quantity 
for  each  station  (usually  313) , 

h  is  the  height  above  the  Earth, 

H  is  the  scale  height,  an  input  quantity  for  each 
station  (usually  7  kilometers), 

n  is  the  index  of  refraction. 

For  the  tropospheric  model,  the  refractive  errors  are  considered 
to  be  independent  of  signal  frequency  since  the  index  of  refraction 
is  virtually  independent  of  frequency  up  to  30,000  MHz, 

In  the  ionospheric  model,  the  index  of  refraction  is  dependent 
upon  more  parameters  than  those  considered  for  the  tropospheric  model. 
The  ionosphere  consists  of  several  belts  of  charged  particles.  The 
F  layer  is  very  much  larger  than  any  other  layer,  and  contains  a 
greater  number  of  charged  particles  than  the  other  layers.  The  F 


(168) 


97 


layer  is  the  one  closest  to  the  Earth's  surfacet  It  Is  subdivided 
into  the  ¥1  and  F2  layers*  In  the  ionospheric  models  the  index 
cf  refraction  is  primarily  dependent  upon  the  height ,  ho#  of  the 
base,  of  the  ionosphere's  F2  layer »  the  maximum  electron  density  of 
the  F2  layer  I  and  the  height  of  the  maximum  electron  density  of  the 
F2  layer. 

Both  index  of  refraction  and  the  height •  hpi  are  dependent  upon 
diurnalp  solar  activity p  seasonal p  and  geographical  variations  as  well 
as  other  miscellaneous  sporadic  variations.  Unlike  the  tropospheric 
model p  the  refractive  errors  in  the  ionospheric  model  are  frequency 
dependent. 

In  constructing  the  model p  the  range  of  the  signal  frequency  has 
been  limited  to  frequencies  above  100  MHZp  since  this  range  of  spec- 
trxim  both  represents  the  situation  of  greatest  Interest  and  enables 
equation  simplification. 

The  relationship  between  the  index  of  refraction  (n)  the  angular 
frequency  of  the  incident  signal  (w) p  and  the  electron  density  in  the 


(169) 


where 

pg  is  the  electron  density  per  cubic  meter p 

e  is  the  electron  charge  (1.6  x  10"^^  Couloumb8)p 

m  is  the  electron  mass  (9.08  x  10"^^  kilograms) 

Co  is  the  permittivity  of  free  space  (8.854  x  10"^^ 
farads/meter) • 

Using  the  first  two  terms  of  the  binomial  expansion  as  an  approx- 
imationp  the  equation  for  the  index  of  refraction  reduces  to: 


n  -  1  -  40.3  ^  10“^^ 

Note  that  f  has  the  dimensions  of  MHz.  This 

2  IT 


(170) 


where  f  - 


equation  is  true  for  frequencies  above  the  critical  frequencyp  f^^p 
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which  is  defined  aet 


fc  -  8.97  po  *  10"*  (MHz) 

where  pq  Is  the  •axlmum  electron  density  per  cubic  meter. 

From  the  definition  of  N  in  Equation  (168),  Equation  (170)  can 
be  written  as 

Pe  -c 
N  -  -  4.03  *  10  * 

The  modal  selected  for  electron  density  versus  height  consists 
of  a  parabolic  variation  below  the  height  of  maximum  electron  density 
matched  to  a  hyperbolic  secant  profile  above  the  maximum.  The  rela¬ 
tionships  are  as  follows! 

Pe  -  Po  [1  -  (1  -  o)^l  for  0  ^  o  ^  ll 

Pe  ■  sach  (o  “  1)]  for  o  ^  1  j 

where 

h  -  ho 
°  “  hn  -  ho 

h  Is  the  height  above  the  Earth,' 

hQ  is  the  height  of  the  base  of  the  F2  layer  (input 
quantity) , 

hiQ  is  the  height  of  the  maxloium  electron  density  In 
the  F2  layer  (Input  quantity) . 

Figure  24  is  a  plot  of  the  ionosphere  model  normalized  with  res¬ 
pect  to  o  and  1/2  (Pg/Po)*  ^n»  PO  parameters  refer  to 

the  ionosphere's  F  layer.  Using  this  model,  the  refractive  effects 
of  the  D  and  E  layers  are  not  singled  out,  because  they  are  quite 
small  in  comparison  with  those  due  to  the  F  layer  and  are  approxi¬ 
mately  accounted  for  by  allowing  the  electron  density  at  the  bottom 
edge  of  the  F  layer  to  be  zero. 


(171) 


(172) 


(173) 
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Figure  24.  Normalized  3-Parameter  Model  of  Atmosphere 


Computation  of  Ray  Bending 

The  major  effect  of  ray  bending  due  to  refraction  Is  that  the 
observed  elevation  angle,  E,  of  a  vehicle  In  view  of  a  station  is 
greater  that  the  actual  or  computed  value  of  elevation,  This 

effect  is  clearly  seen  in  Figure  25*  The  difference  between  the  com¬ 
puted  and  observed  elevation  Is  the  angle  6. 

In  Figure  25  three  vectors  of  Importance  are  shown: 
p,  the  vector  from  the  station  to  the  vehicle, 

p',  a  vector  from  the  station  toward  the  observed 
vehicle  position  (l»e»,  tangent  to  the  ray 
path  at  the  station),  and 

Ap,  a  vector  constructed  perpendicular  to  p* 

The  value  of  azimuth.  A,  and  computed  elevation,  E^,  are  ob¬ 
tained  from  Equations  (156)  and  (157)  using  the  value  of  p,  since 
these  values  do  not  incorporate  refraction  effects.  From  the  geome- 
try  of  Figure  25,  one  can  find  the  components  of  Ap  in  the  east, 
north,  and  vertical  station-oriented  coordinate  system.  Thus: 


(174) 


(175) 


(176) 


Ap  ■  [q  \  tan  6  [T] 


sin  E^  sin  A 
sin  E^  cos  A 
-  cos  Eq 


where  |p|  is  the  magnitude  of  the  vector  p,  and  [T]  is  the  3^3 
matrix  of  Equation  (154)  used  to  transform  the  coordinates  of  Ad 
from  the  station  oriented  coordinate  system  into  the  inertial  system. 
From  the  figure  it  is  also  obvious  that 

p ’  +  Ap  *  p 


or 


p  -  1^1  tan  6[T] 


sin  Ec  sin  A 
sin  E^  cos  A 
-  cos  Eo 


In  the  program  the  vector  p*  is  actually  normalized  to  a  unit  vec¬ 
tor.  The  length  of  p’  is  not  important,  however,  since  it  is  used 
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Figure  25.  Geometry  of  Roy  Potfi  Used  to  Obtain  p 
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only  to  calculate  observed  angle-related  measurements  such  as  azimuth^ 
elevation,  il-cosine,  m-cosine,  X-angle,  Y-angle,  and  indirectly,  hour 
angle  and  declination  as  described  above* 

In  Equation  (176)  all  quantities  are  easily  computed  with  the 
exception  of  6,  the  "error”  in  the  elevation  angle.  The  method  of 
computing  6  used  by  the  program  is  that  due  to  Weisbrod  and  Anderson^ . 
Only  an  outline  of  their  analysis  is  presented  here. 

First,  consider  the  ray  of  Figure  26  entering  an  infinitesimal 
layer  of  thickness  dr  at  an  angle  6*  At  a  boundary  between  two 
infinitesimal  layers  Snell's  law  holds.  Thus; 

n  cos  6  -  constant  (177) 

where  n  is  the  index  of  refraction.  Differentiating  Equation  (177) 
with  respect  to  path  length,  s,  yields  Equation  (178)  below.  Since 
n  varies  only  with  altitude  or  distance  along  an  earth  radius,  r, 
one  obtains: 


^  CO*  S  -  n  4^  sin  6  ■>  0 
as  a  8 

(178) 

dn  dr  .  dQ  ,  „  „ 

—  ^  cos  6  -  n  sin  6  -  0 

(i/‘.0 

but 

^  -  sin  6, 
ds  • 

(180) 

and  from  the  definition  of  curvature  one  gets: 

ds  K 

(181) 

where  <  is  the 

radius  of  curvature  of  the  ray. 

Thus,  after  substi- 

tution.  Equation 

(179)  becones: 

1  1  dn  . 

“  ■>  “  “T"  COS  6  • 

K  n  dr 

(182) 

From  Figure 

26,  an  infinitesimal  length  ds 

of  the  ray  path  is 

given  by: 

dr 

sin  6 

(183) 
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thu8» 


dY 


^  ala  0  * 


Combining  Equations  (18A)  and  (182)  raaults  In  an  axpraaalon  for 


dYI 


dY  ■  “  ^  cot  0  dr 

Also,  from  Figure  26,  It  la  aaaa  that  the  dY'a  of  all  elementary 
layers  are  directly  additive.  Therefore,  the  angle  Yj^  due  to  the 
total  bending  between  layera  bounded  by  the  heights  rj  and  la 

simply . the  Integral  of  Equation  (185) x 


(184) 


(185) 


Yjk 


rk 

i  ^  cot  0  dr 
n  or 


If  a  ray  departs  from  a  given  layer  (where  r  ■  r^^  and  n  ■  n^c)  at  an 
angle  of  0;^,  from  Snell's  law  for  spherical  stratification,  one  hast 


(IbO 


nk  cos  ■  constant 


(187) 


If  the  angle  Yji(  due  to  refractive  bending  Is  computed  for  only 
a  very  small  layer  of  atmosphere  between  the  heights  r  and 


r  -  r^,  the  following  assumptions  are  Justified! 

(1)  ^  "  constant, 

(2)  the  Index  of  refraction  n  Is  close  to  unity,  and 

(3)  Ah  -  rjt  -  rj  «  rj  (Ah  Is  an  infinitesimal  layer 
of  height) • 

Using  these  assumptions  along  with  Equations  (186)  and  (187) ,  Welabrod 
and  Anderson derive  an  expression  for  Yj)^  In  terms  of  the  angle 
0  and  the  refractlvlty  N  (see  equation  (168)  for  Its  definition)^ 
for  the  details  the  reader  Is  referred  to  their  paper.  Their  expres¬ 
sion  for  Yjk  through  a  small  layer  of  atmosphere  Is: 

I'Jk  •  500  6k) 

The  total  bending  through  the  atmosphere  Is  the  sum  of  the 


(188) 
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individual  contributions.  Thus, 


^  N,  1  -  N. 

i"  1  i  ^ 

Y  =  Zu  cnn  rZ - Z - TT - (milliradians) 

500  (tan  +  tan 


(189) 


From  Equation  (187),  the  value  of  cos  P  is  given  as: 

i 


cos  = 


"i-i  ^-1  **1-1 

n.  r . 

1  1 


\-i  ^-1 

n. 

1 


1  - 


Ah 


R  +  h. 
e  1 


(190) 


and  it  follows: 


J  1  "  cos^  p^ 

cos  p.  ~  (191) 

Equations  (189),  (190),  and  (191)  determine  the  value  of  y 
where  P^  is  the  value  of  the  angle  of  incidence  of  the  ray  at  a 
height  h^  =  N(h^)  is  the  refractivity  at  height  h^  obtained 
from  the  atmospheric  model,  and  Ah  is  the  increment  of  height  used 
in  Equations  (189)  and  (190),  Note  that  h^  =  h^  ^  and  that 

h^  is  the  computed  height  of  the  vehicle  obtainable  from  the  program. 

At  the  station  h  =  h^,  which  is  the  height  of  the  station,  and  P^ 

is  set  equal  to  the  computed  elevation,  E^. 

Once  y  has  been  found,  consideration  of  the  geometry  in  Figure 
27  leads  to  the  evaluation  of  6,  the  "error*'  in  the  elevation  angle. 

Thus,  we  note  from  Figure  27: 


e  =  Y  -  (a  -  P)  (192) 

Using  Snell's  law  for  spherical  stratification  and  the  application  of 
the  law  of  sines  to  triangle  SOQ  results  in: 


cos  p 


"o 


nr 


cos  p 


0 


Tq  =  Pq  =  r  cos  a 

Combining  Equations  (193)  and  (19^)  yields: 

R  "O 

P  jj  cos  a  =  cos  (a  -  (a  -  p)) 


(193) 

(19A) 


(195) 
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Figure  27.  Geometry  of  Roy  Path  Used  to  Obtain  3 
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Since  a  -  P  is  a  small  angle,  a  small  angle  approximation  and  an 
approximation  for  is  used^^^^  in  Equation  (195)  to  obtain 

an  expression  for  a  -  P  in  terms  of  refractivity  at  the  station 
(Nq)  and  at  the  vehicle  (N): 

a  -  6  *  (No  -  N)  10“^  cot  P  (196) 

Using  Equations  (192)  and  (196),  the  value  of  c  is  (ietermined • 

Applying  the  law  of  sines  to  triangle  SOP  results  in: 


Combining  Equations 
tan 


cos  (Bo  -  5)  «  r  COL  ((i  +  t )  -  -^ ) 

(194)  and  (197)  gives  an  expression 

5  .  tan  g  +  .0  -  cos  c) _ 

•in  c  +  cos  c  tan  a  -  tan  Bq 


for 


(1^7) 


(197) 


Finally,  using  small  angle  approximations  for  6  and 


c 

tan 

a  + 

£2  i 

MM 

2 

*  "  e  + 

tan 

a  - 

tan  f'o 

(199) 


This  value  of  5  can  now  be  used  in  Equation  (176)  to  find  the 
vector  p',  which  determines  all  angle  related  observations  includin 
the  effects  of  refraction.  Due  to  the  approximations  used  in  the 
analysis,  below  elevations  of  about  5®  the  errors  in  the  propa^^ation 
corrections  amount  to  a  few  percent.  Because  of  this  fact,  propaga¬ 
tion  corrections  due  to  refraction  should  be  computed  only  when  the 
elevation  of  the  vehicle  exceeds  5*. 


Refraction  Effects  on  Range  and  Range  Pxale 

Because  of  signal  retardation  caused  by  the  refractive  grc.di»rt 
of  the  atmosphere,  the  round-trip  range  computation  of  Equation  (158) 
includes  a  correction  term,  Ap.  The  signal  retardation  caused  by  Cxxx 
infinitesimal  layer  of  thickness  dr  is  given  by: 
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■  fi  -  Q5C  g  dr 
^  V  c ' 


-  - 1) 


CSC  S 


dr 


-  N  X  10”^  ^  dr 

C 

The  effect  on  range  between  layers  bounded  by  heights  r  -  rj  and 
r  -  r^  is  given  by: 


{liV) 


r 


Apjk 


’'k 

f 


cdT 


N  X  10“^  CSC  6  dr 

m 

Using  methods  similar  to  the  computation  of  retractive  bending, 
Weisbrod  and  Anderson^ arrive  at  the  expression  for  the  refractive 
round-trip  range  correction,  2ip,  due  to  the  passage  of  the  ray 
through  t\ie  entire  atmosphere: 


(2(1) 


2Ac  -  2  X  10"^ 


y  |Ni_i  +  N^I  (hj^  -  h,_i ) 

sin  6i-i  +  sin 


i-1 


+  10”^ 


^  sin  3i-i  sin  64 

i-nH-i  ^  ^ 


t. 


where 


is  the  refractivity  in  the  i^^  layer, 
is  the  height  of  the  i^^  layer, 

3^  is  obtained  from  equation  (190), 
fj  is  the  transmitted  (or  up)  frequency, 
f2  is  the  received  (or  down)  frequency. 


The  indices  in  the  first  summation  refer  to  layers  in  the  troposphere, 
and  those  in  the  second  refer  to  layers  in  the  ionosphere.  The  values 
of  -  N(hi)  are  obtained  from  the  models  of  refractive  index. 

Due  to  refractive  bending,  the  range-rate  measurement  Of  Equa¬ 
tion  (159)  also  includes  a  correction  term,  Ap,  This  correction 
term  arises  because  the  program  initially  computes  the  component  of 
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vehicle  velocity  along  the  vector  i.e,,  along  the  straight  line 

SP  of  Figure  27*  The  measured  range-rate  is  the  component  of  vehi¬ 
cle  velocity  along  the  tangent  to  the  ray  path  at  the  vehicle,  i.e., 
along  the  straight  line  TP  of  Figure  27. 

The  velocity  of  the  vehicle  and  the  station  is  available  from  the 
program.  If  one  denotes  the  component  of  vehicle  velocity  along  the 
line  SP  by  v^,  and  the  component  of  vehicle  velocity  perpendicular 
to  SP  by  Vy,  then  the  range  rate  computation  without  the  correc¬ 
tion  for  refraction  (i.e.,  Ap  -  0)  is: 


The  measured  component  of  velocity  along  the  tangent  to  the  ray  path 
at  the  vehicle  as  seen  from  Figure  27  is 


pjjj  *  Vx  cos  (y  “  <5)  +  vy  sin  (y  -  6) 
The  correction  to  range-rate,  Ap,  is  given  by: 

Ap  -  Pm  -  Pc 


A  factor  should  also  be  included  to  account  for  the  round-trip  fre¬ 
quency  dependent  effects  in  the  ionosphere.  After  making  a  small  angle 
approximation  the  final  equation  used  for  Ap ,  the  one-way  correc¬ 


tion  to  range-rate,  is: 


Vy  (Y  -  >5) 


where 


Vy  is  the  component  of  vehicle-station  velocity  perpendicular 
to  the  vector  p,  and  lying  in  the  plane  which  includes 
the  center  of  the  Earth, 

f 1  is  the  transmitted  frequency, 

f2  is  the  received  frequency. 


(203) 

(204) 

(205) 

(206) 
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SECTION  III. 


USER'S  GUIDE 

INTRODUCTION 

The  information  in  this  section  is  intended  as  a  guide  for  those 
desiring  to  use  SPACE-A.  SPACE-A  generates  a  vehicle  ephemeris  and 
associated  observations  from  ground  stations  and  has  four  modes  of 
operation.  They  are: 

(1)  Computation  of  vehicle  ephemeris  only; 

(2)  Generation  of  observations  from  ground  stations; 

(3)  Same  as  (2)  with  the  writing  of  the  observations  onto 
tape  in  a  format  suitable  for  processing  by  SPACE-B; 

(4)  Same  as  (2)  with  the  following  included  as  standard 
outputs : 

(a)  the  acquisition  time, 

(b)  total  time  in  sight,  and 

(c)  the  option  of  the  printing  time  of  polar  and 
meridian  crossings. 

DESCRIPTION  OF  INPUT  DATA 

All  input  data  is  read  at  the  beginning  of  each  run.  The  format 
and  arrangement  of  the  data  is  designed  to  make  the  data  prepara¬ 
tion  as  convenient  as  possible.  It  is  seldom  necessary  to  rend  all 
the  data;  in  fact,  most  runs  require  a  minimum  of  cards.  Many  vari¬ 
ables  have  standard  values  which  are  preset  by  the  program  and  should 
be  satisfactory  for  most  cases. 

The  input  data  is  divided  into  three  parts.  The  first  consists 
of  a  single  card,  called  the  basic  input  card,  which  must  be  included 
in  each  run.  The  second  and  third  parts  are  designated  group  1  and 
group  II  inputs,  respectively.  Each  group  is  divided  into  several 
sections,  each  of  which  must  be  preceded  by  a  card  containing  the 
section  number,  in  integer  format,  in  columns  four  and  five.  Any 
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section  may  be  omitted  from  the  input  deck*  Parts  of  sections,  how¬ 
ever,  must  not  be  omitted*  See  Figure  28  for  the  setup  of  a  sample 
input  deck* 

The  basic  input  card,  which  precedes  each  case,  contains  four 
variables*  The  first  pertains  to  the  standard  values  which  may  be 
set  in  place  of  the  group  II  data  and  will  be  explained  later*  The 
second,  MDE,  indicates  the  function  to  be  performed  (i*e*,  trajectory 
computation,  observables,  etc*)*  The  third  indicates  the  desired  pre¬ 
cision  and  determines  the  set  of  standard  values  used  for  integration* 
The  last  specified  Encke  or  Cowell  integration* 

The  group  I  inputs  consists  of  several  sections*  Any  of  these 
sections  may  be  omitted  from  the  input  deck  if  the  variables  contained 
in  the  section  are  not  needed  for  the  run  or  in  the  case  of  stacked 
cases,  the  values  from  the  last  case  are  to  be  used  again*  The  des¬ 
cription  of  each  variable  given  in  the  data  summary  (see  Table  IV)  is 
sufficient  to  understand  most  of  the  inputs*  Therefore,  only  a  brief 
discussion  of  the  variables  given  in  Table  IV  is  presented  here. 

In  the  group  I  input  data,  if  the  variable  KLM2  of  Section  2  is 
0  or  1,  the  initial  conditions  are  assumed  to  be  referred  to  the  base 
date  system  or  true  system  of  date,  respectively* 

Sections  7  and  9  deal  with  observations  from  observing  systems 
on  vehicles  and  with  powered  flight  and  have  been  omitted  since  they 
are  not  presently  used*  The  variables  in  Section  8  control  the  out¬ 
put  and  are  discussed  in  DESCRIPTION  OF  OUTPUT  INFORMATION* 

The  end  of  the  group  I  input  data  is  indicated  by  a  card  con¬ 
taining  a  10  in  columns  four  and  five*  This  card  must  always  be  in¬ 

cluded  in  the  data  deck*  After  reading  this  card  the  program  advances 
to  the  group  II  inputs* 

Before  the  program  reads  the  group  II  input  data,  the  val»e  of 
the  first  variable  contained  on  the  basic  input  card,  KSTDRD,  is 

tested*  If  this  variable  is  negative,  standard  values  of  the  vari- 
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Figure  28.  Sample  Input  Deck 


ablts  in  the  group  II  input  data  are  computed  from  stored  information 
and  the  program  then  begins  to  read  the  group  II  inputs  to  augment  or 
replace  any  of  the  standard  values*  Table  V  gives  a  list  of  the  stan¬ 
dard  values  used  by  the  program.  If  KSTDRD  is  not  negative ^  none  of 
the  standard  values  is  set  and  all  values  should  be  read  in  unless 
values  from  a  preceding  stacked  case  are  used.  The  end  of  group  II 
inputs  is  indicated  by  reading  a  card  with  a  20  in  integer  format  in 
columns  four  and  five. 


114 


Table  IV. 

Summary  of  Input  Data 
Basic  Input  Card 


COLUMN  NAME  FORMAT 


DESCRIPTION 


2-5  KSTDRD  15 

6-10  .  MDE  15 


11-15  PRECIS  F5.0 


16-20  CEPID  F5.0 


<  0  Will  want  at  least  some 
standard  Inputs  from 
group  II 

^  0  No  standard.  All  values 

will  be  read  in. 

■  1  Trajectory  computation. 

•  2  Observable  computations 

without  simulated  data. 

*  3  Observable  computations 

with  simulated  data. 

■  A  Prediction  mode 

-  1  Low  precision  level. 

■  2  Intermediate  precision 

level. 

"  3  High  precision  level 

■  0,1  Encke  integration 
■>  -  1  Cowell  integration 
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Group  I  Inputs 


SECTION  CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

1  1 

2-72 

ITITLE(1-12) 

12A6 

Title 

2 

2-5 

NYEARP 

15 

Year  of  Initial 
conditions 

6-10 

DAYS 

F5.0 

Day  of  Initial 
conditions. 

11-15 

HRS 

F5.0 

Hour  of  Initial 
day. 

16-20 

HMIN 

F5.0 

Minute  of  Initial 
hour. 

21tA0 

SEC 

E20.16 

Seconds  of  Initial 
minute. 

3 

1-2A 

TMAX 

E2A.16 

Time  of  run  In 
hours. 

2  1 

2-5 

KLM 

15 

Indicator  for  units 
of  position  and 
velocity  vector. 

-  1  KM, KM/SEC 

-  2  ER,ER/HR 

-  3  FT,Fr/SEC 

-  A  MI, MI /HR 

-  5  NM,NM/HR 

-  6  NM, FT/SEC 

-  7  AU,AU/HR 

6-10 

KLMl 

15 

Indicator  for  coor¬ 
dinate  system  of 

Input  vectors. 

■  1  Cartesian  coord. 

WE  -  0 

■  2  Cartesian  coord. 

Compute  WE 

-  3  Geodetic  long, 

lat,  alt, 
Ve.Vn.Vh;  WE  -  0 

-  A  Geodetic  long, 

lat,  alt, 

Ve.Vn.Vh; 

Compute  WE 
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SECTION  CARD 


Grou£iJ^J[ngut£ 

(continued) 

COLUMNS  NAME  FORMAT 


11-15  KLM2  15 


16-20  KSNAP  15 


21-25  KLM3  15 


DESCRIPTION 

■■  5  Geodetic  long, 
lat,  alt, 

|v| ,  Y,  az; 

WE  -  0 

■  6  Geodetic  long, 

latp  alt^ 

|v| ,  Y»  az; 
compute  WE 

■  7  Geocentric  ra, 

declf  altt 

^ra  t^d  f t 

WE  -  0 

■  8  Geocentric  ra, 

declp  altp 
VratVdfVh; 
compute  WE 

■  9  Geocentric  ra, 

decl^  altp 
|v| ,  6,  az; 

WE  -  0 

«  10  geocentric  ra, 
decl,  alt, 

|v| ,  6,  az; 
compute  WE 

Indicator  for  nuta¬ 
tion  and  precession 
of  input  vector, 

-  0  no 

-  1  yes 

Indicator  for  nuta¬ 
tion  and  precession 
of  vectors  during 
run, 

■  0  no 

■  1  yes 

Indicator  for  libra- 
tion  of  input  vectors, 

-  0  no 

-  1  yes 
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Group  I  Inputs 

(continued) 

SECTION 

CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

26-30 

KLIBR 

15 

Indicator  for  libra- 
tion  of  vectors 
during  run 

■  0  no 

■  1  yes 

31-35 

MRREF 

15 

Indicator  for  ini¬ 
tial  reference  body 

-  1  Earth 

■  2  Sun 

■  3  MOon 

-  A  Mars 

■  5  Venus 

■  6  Jupiter 

-  7  Saturn 

2 

1-24 

RCIN(l) 

E24.16 

Initial  position 
vector 

25-48 

(2) 

E24,16 

See  KLM  and  KLMl 
for  units  and  type* 

(3) 

E24.16 

3 

1-24 

RCIN(l) 

E24.16 

Initial  velocity 
vector. 

25-48 

(2) 

E24.16 

See  KLM  and  KLMl 
for  units  and  type. 

49-72 

(3) 

E24.16 

3 

1 

2-5 

PASFX 

F5.0 

Total  number  of 
passes. 

4 

1 

2-5 

KS2BY 

15 

Indicator  for  two- 
body  integration 

■  0  no 
-  1  yes 

6-10 

KSPLT 

15 

Indicator  for  in¬ 
clusion  of  planetary 
perturbations 

-  0  no 
•  1  yes 
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SECTION 


Group  I  Inputs 
(continued) 

CARD  COLUMNS  NAME  FORMAT 

11-15  KSOBL  15 


16-20  KSDRG  15 


21-25  KSRAP  15 


26-30  KSDRGM  15 


31-35  KSDRGV  15 


36-40  KSMNOB  15 


41-45  KRF  15 


DESCRIPTION 

Indicator  for  inclu¬ 
sion  of  oblateness 
perturbatlcn# 

■  0  no 

■  1  yes 

Indicator  for  inclu¬ 
sion  of  Earth  drag 
perturbation* 

■  0  no 

■  1  yes 

Indicator  for  inclu¬ 
sion  of  radiation 
pressure  perturba¬ 
tion* 

■  0  no 

■  1  yes 

Indicator  for  inclu¬ 
sion  of  Mars  drag 
perturbation* 

■  0  no 

-  1  yes 

Indicator  for  inclu¬ 
sion  of  Venus  drag 
perturbation* 

■  0  no 

-  1  yes 

Ind i ca tor  f or  inclu¬ 
sion  of  Moon  oblate¬ 
ness  perturbation* 

■  0  no 

■  1  yes 

Indicator  for  inclu¬ 
sion  of  refraction 
effects* 

■  0  no 

■  1  yes 
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Group  I  Inputs 


(continued) 

SECTION  CARD  COLUMNS  NAME  FORMAT  DESCRIPTION 

46-50  KECLPS  15  Indicator  for  inclu¬ 

sion  of  eclipse  in¬ 
formation  print. 


-  0  no 

-  1  yes 


5  1 

2-5 

MAXSTA 

15 

Total  number  of 
stations  used  in 

run. 

2 

2-5 

K 

15 

Station  number. 

10-15 

STANM(L) 

A6 

Station  name. 

20-30 

TYPE(L) 

4X111 

A(1)+1,E2*A(2) 

+l,E4*A(3)4*i,E6*A(4) 
+1,E8*K  where 
A  ■  observation 
types  used  by 
station  K  in 
ascending  order. 
L  ■  packed  station 
number 


3 

1-24 

STALT(K) 

E24.16 

Latitude  of  station 

K, 

25-48 

STALN (K) 

E24.16 

Longitude  of  station 

K. 

49-72 

START (K) 

E24.16 

Altitude  of  station 

K. 

4  , 

1-36 

RRATE(1,L) 

E36.16'^ 

r Repetition  rates 

1  (hrs). 

37-72 

RRATE(2,L) 

E36.16 

\  \  For  each  observation 

5 

1-36 

RRATE(3,L) 

E36.16 

Llype 

37-72 

RRATE(4,L) 

E36.16J 

6 

1-18 

TDELAY(1,L) 

E18.8 

r Times  in  hrs  before 

19-36 

TDELAY(2,L) 

E18.8 

J  which  each  observa- 

37-54 

TDELAY(3,L) 

E18.8 

\  tion  is  not  to  be 

55-72 

tdelay^.l) 

E18.8 

V  computed. 
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Group  I  Inputs 
(continued) 


SECTION  CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

7 

1-18 

FUP(K) 

E18.8 

Station  transmit 
freq,  (MHz), 

19-36 

FDOWN(K) 

E18,8 

Station  receive 
freq,  (MIIz) , 

8 

1-2A 

STA0R(NCDST+1) 

E24.16 

AEE  station 
rotation  angle. 

25-48 

STA0R(NCDST+2) 

E24.16 

AEV  station 
rotation  angle. 

49-72 

STAOR(NCDST+3) 

E24.16 

AEN  station 
rotation  angle. 

9 

1-24 

STA0R(NCDST+4) 

E24.16 

4U 

^station  lo¬ 
cation 

25-48 

STA0R(NCDST+5) 

E24.16 

.V 

errors 
caused  by 

49-72 

STAOR(NCDST+6) 

E24.16 

AW 

geodetic 
^net  error. 

10 

1-24 

STA0R(NCDST+7) 

E24.16 

NO;  refractlv- 
vlty  at  sea 
level. 

25-48 

STA0R(NCDST+8) 

E24.16 

H;  troposphere 
scale  fact. 

49-72 

STA0R(NCDST+9) 

E24.16 

PO;  max,  elec¬ 
tron  density. 

11 

1-24 

STAOR(NCDST+10) 

E24.16 

HO; 

of 

loser  limit 
ionosphere. 

25-48 

STA0R(NCDST+11) 

E24.16 

HM;  Ht  of  PO 
(KM) 

49-72 

STAOR(NCDST+12) 

E24.16 

-  open  - 

12 

1-24 

STA0R(NCDST+13) 

E24.16 

AT 

for  timing. 

25-48 

STA0R(NCDST+14) 

E24.16 

Bias  added  for 
obser ,A(1) 

49-72 

STA0R(NCDST+15) 

E24.16 

Bias  added  for 
obser ,A(2) , 

121 


Group  I  Inputs 
(continued) 


SECTION  CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

13 

1-24 

STAOR(NCDST+16) 

E24.16 

Bias  added  for 
obser • A(3) 

25-48 

STAOR(NCDST+17) 

E24.16 

Bias  added  for 
obser, A(4) 

REPEAT 

CARDS  2-13  FOR  EACH 

STATION 

6  1 

1-18 

DAREAS 

E18.8 

Effective  sur¬ 
face  area  of 
vehicle  per¬ 
taining  to 
drag 

19-36 

PAREAS 

E18.8 

Effective  sur¬ 
face  area  per¬ 
taining  to 
radiation  pres¬ 
sure  (ft^). 

37-54 

SPADD(6) 

E18.8 

Mass  of  vehicle 
(Ib-tnasses) 

7 

OPEN 

8  1 

2-5 

lUNIT 

15 

Indicator  for 
output  units 
(see  KLM) 

6-55 

IPSEC(I),I-1.10 

1015 

Indicator  for 
suppression  of 
each  of  10  print 
sections. 

2 

1-18 

19-36 

DTP  I 

DTSUP 

E18.8 

E18.8 

Print  portion 
(hrs)  and  sup¬ 
press  portion 
(hrs),  of  total 
print  period. 

37-54 

PRATE 

E18.8 

Interval  with¬ 
in  DTPI  at  which 
to  print. 

9  OPEN 

10  END  OF  GROUP  I  INPUTS 
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Groug_^I_^tJ2U5£ 


SECTION  CARD 

COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

1  1-7 

1-72 

(DT3(I»J), 
1-1,3)  J-1,7 

3E2A.12 

Integration  in¬ 
tervals  for  each 
of  7  working 
bodies  for  near, 
medium  and  far 
reference. 

2  1 

1-72 

Rl(l-3) 

3E2A,12 

R1  and  R2  are 

2 

1-72 

Rl(A-6) 

3E2A.12 

distances  in 

E.R,  for  each  of 

3 

1-72 

Rl(7) ,R2(l-2) 

3E2A.12 

7  working  bodies 

A 

1-72 

R2(3-5) 

3E2A.12 

for  switching 
from  near  to 

5 

1-A8 

R2(6-7) 

3E2A.12 

medium  and  me¬ 
dium  to  far  in¬ 
tegration  inter¬ 
vals* 

3  1 

1-2A 

RTl 

E2A.12 

Values  used  as 

25-A8 

RT2 

E2A.12 

tolerances  in 

rectification 

criteria. 

A  1 

1-18 

DHl 

E18.8 

Troposphere 
integration 
step  size  (KM)* 

19-36 

DH2 

E18.8 

Ionosphere  in¬ 
tegration  step 
size  (KM)* 

37-5A 

H2 

E18.8 

Upper  limit 
troposphere 
(KM). 

55-72 

HA 

E18.8 

Upper  limit 
ionosphere  (KM) . 

5  1 

2-5 

KOBLAT 

15 

Number  of  ob¬ 
lateness  coeffi¬ 
cient  terms* 

2 

2-5 

N 

15 

N  index* 

6-10 

M 

15 

M  index* 

3 

1-2A 

SORCl 

E2A.12 

Value  of  C 
coefficient* 
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Grou£^^IIn£ut^ 


(continued) 

SECTION  CARD  COLUMNS 

NAME 

FORMAT 

DESCRIPTION 

25-48 

SORC2 

E24.12 

Value  of  S 
coefficient 

REPEAT  CARDS  2-3  UNTIL  KOBL^T  VALUES  HAVE  BEEN  READ  IN. 


6 

1 

2-5 

MBMAX 

15 

Number  of  working 
bodies. 

6-10 

KWBMUd) 
1*1, MBMAX 

15 

Indices  of 
working  bodies. 

2 

1-72 

TPMAT8(1) 
I-l, MBMAX 

3E24.12 

Gravitational 
constants  of 
working  bodies. 

REPEAT 

CARD  2 

FOR  EACH  VALUE 

OVER  3N  NEEDED. 

7 

1 

1-24 

DYN(48) 

E24.12 

Solar  flux  in 
10-22  w/m^-Hz 
at  10.7  cm. 

25-48 

DYN(49) 

E24.12 

Open 

8 

1 

1-72 

DYN(31-53) 

3E24.12 

Coefficients 
for  lunar  ob¬ 
lateness 
(kg  •  km^). 

9 

1 

1-24 

COMB(l) 

E24.12 

Velocity  of 
light. 

10 

1 

1-72 

PRNT3(l-3) 

3E24.12 

Print  intervals 
(hrs)  for  near 
medium  and  far 
reference. 

11 

1 

1-24 

EMIN 

E24.12 

Minimum  value 
of  elevation 
angle  (radians) . 

12 

1 

1-18 

RTO 

E18.8 

Ratio  of 

Nordsieck  inte¬ 
gration  interval 
to  that  in  Runge- 
Kutta. 
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SECTION 

13 


14 


15 


16 

17 

18 

19 

20 


Group  II  Inputs 
(continued) 


CARD 

COLUMNS 

NAME 

1 

1-24 

DSPL 

1 

1-72 

RATEC( 1-3,1) 

2 

1-72 

RATEC(l-3,2) 

1-10 

1-72 

XMACH(I) 

1-1,40 

11-20 

1-72 

CDT(I) 

1-1,40 

FORMAT 

DESCRIPTION 

E24.12 

Special  inte¬ 
gration  inter¬ 
val  in  A4  mode 
to  obtain  ac¬ 
quisition  time 
(hrs). 

3E24.12 

Rotation  vector 
used  in  Mars 
drag  computa¬ 
tions. 

3E24,12 

Rotation  vector 
used  in  Venus 
drag  computa¬ 
tions. 

4E18.8 

Mach  number 
table. 

4E18.8 

Drag  coefficient 
table 

OPEN 

OPEN 

OPEN 

OPEN 

END  OF  GROUP  II  INPUTS 
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Table  V. 


Standard  Values  for  Group  II  Inputs 


Section  1. 

DT3  ARRAY  (3x7 

DT3(1,I);  1-1,3, A, 5, 6, 7 
DT3(2,I);  1-1,3, A, 5, 6, 7 
DT3(3,I);  1-1,3, A, 5, 6, 7 
DT3(1,2)  -  A. 

DT3(2,2)  -  6. 

DT3(3,2)  -  10. 


array  of  integration  Intervals) 
.1953125  X  10-2 
.15625  X  10"! 

.125 


Note:  The  above  values  are  always  used  for  Encke  integration  and  for 
Cowell  when  low  precision  is  specified.  When  Cowell  and  inter¬ 
mediate  precision  is  specified,  each  value  in  the  DT3  array  is 
divided  by  2;  when  Cowell  and  high  precision  is  specified  each 
value  is  divided  by  A.  The  precision  desired  is  controlled  by 
the  variable  CEPID  on  the  Basic  Input  Card. 


Section  2. 

R1  ARRAY  (1  X  7  array  of  reference  body  change 

criteria) 

R1(I)  -  A  •  RADII(I) 

R2(I)  -  10  •  RADII (I) 

where  RADII (I)  is  the  radius  of  the  Ith  working  body  (Earth,  Sun, 
Moon,  Venus,  Mars,  Jupiter,  Saturn). 

Section  3. 

RTl  -  1.  X  10”^ 

RT2  -  1.  X  10"^ 

If  low  precision  is  specified 
RTl  -  1  X  10"** 

RT2  -  1  X  10-** 
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Table  V. 

(continued) 


Section  4. 

DHl  -  1. 

DH2  -  10. 

H2  -  40. 

H4  -  2,000. 

Section  5. 


KOBLAT  -  3 


C20  ■  - 

1082.3  * 

10”® 

S20 

C30  ■ 

X 

• 

CM 

10”® 

S30 

C40  ■ 

1.8  * 

10”® 

S40 

Section  6. 

KBMAX  -  7 

KWBMU(I)  -  I;  I  -  1,2, 3, 4, 5, 6, 7 


TPMAT8(1) 

- 

19.909416 

(2) 

- 

6629965.8 

(3) 

- 

.  .24478289 

(4) 

- 

16.1983009 

(5) 

- 

2.14364682 

(6) 

- 

6338.16258 

(7) 

- 

1897.36734 

Section  7. 

DYN(48)  -  200. 

Section  8. 

DYN(51)  -  .88746  *  kg  •  ko^ 

DYN(52)  -  .88764  «  kg  •  k«2 

DYN(53)  -  .88807  *  10^^  kg  *  ka^ 
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Table  V. 
(continued) 


Section  9. 

COMB(l)  -  169210, 5800A 

Section  10. 

PRNT3(1)  -  .25 
PRNT3(2)  -  .5 
PRNT3(3)  -  1. 

Section  11. 

EMIN  -  (5.  degrees;  expressed  In  radians) 

Section  12. 

RTO  -  3 


Section  13. 

DSPL  -  1.  X  10"3 


Section  14. 

RATEV(1,1) 

RATEV(2,1) 

RATEV(3,1) 

RATEV(1,2) 

RATEV(2,2) 

RATEV(3,2) 


0.0 


0.0 

0.255175469 

0.0 

0.0 


0.0 
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Table  V. 
(concluded) 


Section  15. 

MACH  TABLE  ARRAY  (XMACH(I) ,  I  -  1,40) 


.5 

1.5 

25. 

96. 

0.0 

1.85 

30. 

97. 

.5 

2.0 

50. 

98. 

.75 

■2.7 

70. 

99. 

.85 

3.4 

90. 

100. 

.9 

4.2 

91. 

101. 

1.04 

5.6 

92. 

102. 

1.1 

6.75 

93. 

103. 

1.2 

8.5 

94. 

104. 

1.3 

15. 

95. 

105. 

DRAG  COEFFICIENT  TABLE  (CDT(I) 

.8 

1.4 

1.14 

1.14 

.8 

1.38 

1.14 

1.14 

.82 

1.36 

1.14 

1.14 

.92 

1.285 

1.14 

1.14 

.99 

1.23 

1.14 

1.14 

1.04 

1.19 

1.14 

1.14 

1.18 

1.16 

1.14 

1.14 

1.21 

1.155  ■ 

1.14 

1.14 

1.26 

1.15 

1.14 

1.14 

1.3 

1.45 

1.14 

1.14 

40) 
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DESCRIPTION  OF  OUTPUT  INFORMATION 

The  input  quantities  which  determine  the  form  of  the  output  of 
the  program  are  given  in  Section  8  of  the  Group  I  inputs  in  Table  IV. 
The  choice  of  output  units  are  determined  by  lUNIT,  The  choices  of 
units  with  corresponding  values  of  lUNIT  are  given  below. 

Value  of  lUNIT  OUTPUT  UNITS 

1  KM, KM/SEC 

2  ER,ER/HR 

3  FT, FT/SEC 

4  MI, MI /HR 

5  NM,NM/HR 

6  MM, FT/SEC 

7  AU,AU/HR 

Time  may  be  divided  into  periods  of  print  and  suppression  of 
print  by  setting  the  input  variables  DTPI  and  DTSUP.  A  period  of 
printing  (DTPI)  is  executed  followed  by  a  period  of  suppression  (DTSUP) 
and  the  cycle  is  then  repeated. 

The  variable  PRATE,  is  used  to  determine  the  print  interval  with¬ 
in  DTPI.  Since  conditions  occur  where  the  print  routine  is  entered 
more  than  once  with  the  same  time,  printing  is  normally  only  executed 
upon  the  first  entry.  However,  by  using  a  negative  PRATE,  printing 
will  occur  each  time  the  print  subroutine  is  entered  during  the  print 
interval  determined  by  DTPI. 

The  user  has  the  option  of  choosing  the  information  he  would  like 
printed  by  setting  IPSEC(I)  »  1,  where  I  is  the  output  section  number 
given  in  Table  VI,  that  contains  certain  information.  By  reading 
a  zero  for  IPSEC(I),  the  printing  of  the  I^h  section  is  surpressed. 
There  are  ten  sections  of  output  information  as  shown  in  Table  VI.  See 
Figure  29  for  a  sample  listing  with  various  sections  labeled. 
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Table  VI. 

Ten  Sections  of  Output  Information 


Section  1. 

Program  time  in  hours  since  launch; 

Current  program  date  in  hours,  minutes,  seconds, 
day  and  year; 

Planetary  body  which  is  the  current  reference. 

Section  2. 

Station  number; 

Observation  types  which  station  handles. 

Section  3. 

Components  and  magnitudes  of  position  and  velocity 
vectors  in  units  specified  by  the  user. 

Section  4. 

Components  and  magnitudes  of  perturbations  in 
position,  velocity  and  acceleration  in  units 
specified  by  the  user. 

Section  5. 

Components  and  magnitudes  of  vectors  between 
vehicle  and  each  planetary  body  in  specified 
units. 

Section  6. 

Components  and  magnitudes  of  vectors  between 
the  Earth  and  each  of  the  other  planetary  bodies 
in  specified  units. 

Section  7. 

Components  and  magnitudes  of  vectors  between 
the  Sun  and  each  of  the  planetary  bodies 
in  specified  units. 
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Table  VI. 
(concluded) 


Section  8. 

Right  ascension; 

Declination; 

Greenwich  hour  angle; 

Longitude  and  latitude  of  Earth  subsatelllte  point 
Geocentric  altitude^  azimuth  and  elevation; 
Geodetic  azimuth  and  elevation* 

Section  9* 

True^  eccentric  and  mean  anomalies; 

Semi-major  axis; 

Eccentricity; 

Inclination; 

Argument  of  perigee; 

Perigee  passage  time; 

Right  ascension  of  ascending  node; 

Period; 

Mean  motion; 

Perigee  and  apogee  heights; 

Unit  vector  to  perigee; 

Unit  angular  momentum  vector* 

Section  10* 

Station  name; 

Program  time  In  hours; 

Station  location; 

Observation  values* 
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PROGRAM  DECK  SETUP 

This  section  describes  the  deck  setup  for  running  the  A-mode  of 
SPACE  on  the  7030  computer#  The  program  is  executed  under  the  control 
of  MCP  (Master  Control  Program)#  The  binary  cards  containing  the  com¬ 
piled  routines  are  kept  on  a  magnetic  tape  at  the  MITRE  7030  computer 
facility  and  only  a  few  control  cards  and  the  data  deck  are  needed  by 
the  user# 

Those  cards  shown  in  Figure  30  that  have  a  B  or  T  punched  in 
column  1  are  control  cards  for  the  MCP.  They  are  described  in  greater 
detail  in  the  MITRE  7030  Facility  Manual. 

JOB  Card 

The  job  card  is  always  the  first  card  of  an  input  deck.  This 
card  contains  necessary  information  for  accounting  and  no  job  can  be 
run  without  it#  The  fields  on  the  JOB  card  are  variable  in  length 
and  separated  by  commas# 


Format 
1  10 

B  JOB, 'JOBID' , 'NAME' , ’PROJECT’ DEPT*, ’MAXTIME’ ,'DEST* 


Field  Name 

Description 

B 

Punched  in  column  1 

JOB, 

Punched  in  columns  10-13 

’JOBID' 

1  to  7  character  job  identic 
f ication 

'NAME' 

Name  of  programmer  (up  to  13 
characters) 

' PROJECT ' 

Project  number  (3  to  9 
characters) 

'DEPT' 

Department  number  (3 
characters) 

'MAXTIME' 

Maximum  running  time  in 
minutaa  (up  to  3  charactera) 

'DEST' 

Destination  of  output 
(up  to  6  characters) 
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Example: 

1  10 

B  JOB ^ SPACE, BOND  J,  007,  D8A,  13,  AlAl 

TYPE  Card 

This  card  defines  the  type  of  job  (GO),  and  the  proper  monitor 
system  (FORTRAN)  to  the  MCP  to  execute  this  program. 


Format: 

1  10 

B  TYPE, GO, FORTRAN 

ADDIO  Card 

This  card  contains  the  library  number  of  the  program  tape  that 
contains  the  binary  decks  of  the  SPACE  program.  Since  modifications 
to  the  program  may  be  made  in  the  future,  the  user  should  periodically 
check  the  current  tape  number  with  Department  D-8A. 


Format: 

10 

ADDIO, PLBXXXX 


1 

BLIBl 


Field  Name 
BLIBl 

ADDIO, PLBXXXX 


SUBTYPE, BIN  Card 

This  card  tells  the  MCP  that 


Description 

Punched  in  columns  1-5 

Beginning  in  column  10, 
the  XXXX  indicate  the 
field  where  the  tape 
number  is  to  be  punched. 


a  cards  that  follow  are  binary. 


Format : 

1  10 

T  SUBTYPE, BIN 

TMAIN  Card 

This  card  results  in  the  construction  of  a  dummy  main  program 
which  calls  the  executive  routine  of  SPACE  (EXECA).  EXECA  is  compiled 
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aa  a  subroutine  and  the  binary  cards  are  contained  on  the  library  tape« 
This  card  Is  the  equivalent  to  the  Fortran  routine* 

CALL  EXECA 

RETURN 
END 


Format: 

1  10 

T  MAIN, EXECA 

Processed  FIOD  Deck 

Since  the  programmer  cannot  choose  the  Input-output  channels  to 
be  used  (they  are  assigned  by  MCP  at  execution  time),  the  FIOD  deck 
specifies  the  I/O  channel  requirements  of  the  program*  The  processed 
FIOD  deck  Is  the  result  of  compiling  the  following  subprogram: 


1 

10 

T 

SUBTYPE, FIOD 

B2 

IOD,$READER 

B3 

IOD,$PRINTER 

B9 

lOD.TAPE,,,,,, 

,SAVE 

B8 

lOD.TAPE . . 

,SAVE 

B 

REEL,PLBDL680 

END 

Logical  units  2  and  3  are  the  system  Input  and  output  units  res¬ 
pectively.  Logical  unit  9  Is  used  for  the  output  of  the  observation 
tape  generated  when  that  option  Is  exercised.  Logical  unit  8  Is  used 
for  the  Input  of  the  ephemerls  tape  containing  the  positions  of  the 
Moon  and  planets*  The  reel  card  specifies  that  the  ephemerls  tape, 
whose  label  Is  PLBL680,  Is  to  be  mounted* 

SUBTYPE. DATA  Card 

This  card  tells  the  MCP  that  the  following  cards  contain  the 
data  for  the  program* 
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Format! 


1  10 

T  SUBTYPE, BIN 

Input  Data  Deck 

This  Is  the  data  deck  described  in  the  Description  of  Input  Data. 


Figure  30.  Deck  Setup 


139 


SECTION  IV. 
PROGRAMMING 


GENERAL  FLOW  OF  PROGRAM 

The  general  flow  of  the  program  is  shown  in  Figure  31*  The  user 
does  not  need  to  specify  the  number  of  cases  to  be  run^  as  the  program 
returns  to  read  the  data  for  the  next  case  until  all  data  cards  are 
exhausted. 

The  program  structure  and  linkage  of  the  program  between  basic 
routines  is  shown  in  Figure  32.  It  is  to  be  noted^  however^  that  the 
figure  does  not  indicate  every  linkage  but  only  the  major  ones.  For 
instance,  in  some  cases  it  is  necessary  for  05SERA  (which  usually  only 
computes  observations)  to  call  the  integration  package. 

The  routines  listed  under  the  heading  of  General  Routines,  in 
Figure  32  are  used  by  many  of  the  major  routines.  The  linkage  of 
these  routines  is  not  included  as  this  would  defeat  the  purpose  of 
the  figure,  which  is  to  indicate  the  position  of  the  major  routines 
in  the  overall  structure  of  the  program. 

DESCRIPTION  OF  SUBROUTINES 

The  following  part  of  the  document  gives  a  brief  description  of 
the  subroutines  used  by  SPACE-A.  Each  subroutine  description  gives 
the  purpose  and  a  brief  description  of  the  subroutine,  as  well  as  the 
other  subroutines  which  call  and  are  called  by  the  given  subroutine. 
For  understanding  the  details  of  each  subroutine  one  ought  to  refer 
to  the  pertinent  portions  of  Section  II  and  Section  III  dealing  with 
the  function  of  that  subroutine  as  well  as  referring  to  the  computer 
listing. 

It  should  be  noted  that,  as  of  the  writing  of  this  document, 
certain  subroutines  have  not  been  tested.  These  non-operational  sub¬ 
routines  ares  LUNOBL,  MVDRAG,  OBD,  PFINIT,  PFLGHT,  STACUL.  These 
subroutines  deal  with  lunar  oblateness  gravitational  perturbations. 
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BEGIN 


RETURN  TO 
FIND  NEXT 
ACTIVITY 
TIME 


RETURN 
TO  READ 
DATA  FOR 
NEXT  CASE 


Figure  31 .  General  Flow  of  Program 
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1-872 


EXECA 

(EXECUTIVE  ROUtiNE) 


INPUT  A 

(READS  DATA) 


MAIN  A 

(CONTROLS  MAIN 
FLOW  OF  PROGRAM) 


XFORM 

(READS  a  CONVERTS 
INITIAL  CONDITIONS) 


TIMNGA 

(FIND  NEXT 
ACTIVITY  TIM^ 


OBSER  A 

(COMPUTES 
OBSERVATIONS) 


PFINIT 

(INITIALIZE 
POWER  ED  FLIGHT) 


STACUL 

(CHECK  FOR 
OCCULATION) 


PRINT  A 

(PRINT  SELECTED 
INFORMATION) 


COWELL 


INTEGRATION 
ROUTINES 


OBD  (ON-BOARD  OBSERVATIONS) 


•  ATIM  (MODE  4  COMPUTATION) 


STAPOS  (STATION  VECTOR) 


ENCKE 


♦  MODELA  (INDEX  OF  REFRACTION) 


CITGRA  ( CONTROLS  INTEGRATION) 
CCHREF  (  CHANGE  REFERENCE  BODY) 
CRSTRE  (  SAVE /RESTORE  VECTORS) 
CINT  (INTEGRATION) 


Bl"^)' 


GENERAL  ROUTINES 

DDOT  (COMPUTE  A- B^) 

DMTML  (COMPUTE  [aITbI/IAI 

DOMUD  (CHECK  FOR  ERROR) 

EPHEM  (READ  EPHEMERIS  TAP^i 

FIX  (UNPACK  WORD) 

NUTPRE  (COMPUTE  NUTAT./ 
PREC.) 

SERVCE  (COMPUTE  AX  B, 

ICI ) 


EITGRA  (CONTROLS  INTEGRATION) 
ECHREF  (  CHANGE  REFERENCE  BODY) 
ERSTRE  (SAVE/RESTORE  VECTORS) 
EINT  (INTEGRATION) 

KEPLER  (NOMINAL/2-BODY  ORBIT) 
RECT  (RECTIFICATION  OF  ORBIT) 


I 


DERIV 

(COMPUTES  ACCELERATIONS) 


I 


OBLDRG 

(OBLATENESS 
AND  DRAG) 


LUNOBL 

(LUNAR 
OBLATE NESS) 


MV DRAG 

(MARS/VENUS; 


DENSTY 

(COMPUTES 
DENSITY 
OF  AIR) 


ECLIPSE 

PFLGHT 

(CHECK  a  PRINT 

(POWERED 

ECL.  INFO.) 

FLIGHT) 

Figure  32.  Program  Structure 
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Mars  and  Venus  drag  model p  on-board  obiervattona ^  parturbatlona  due 
to  vehicle  thruatp  and  observatlona  dealing  with  occultation. 

EXECA 

Purpose; 

This  is  the  executive  routine  for  SPACE~A« 

Description; 

The  program  alternates  between  calling  INPUTA  and  MAINA  until 
all  cases  have  been  exhausted. 

Subroutines  Called:  INPUTA,  MAINA 

ATIM 

Purpose: 

This  subroutine  is  used  in  the  visibility  computation  mode  (mode 
4).  It  determines  the  time  a  vehicle  comes  within  sight  of  a  given 
ground  station  or  the  time  of  polar  baseline  crossing  (the  time  that 
the  vehicle  crosses  the  east-west  vertical  plane  through  the  zenith) 
and  meridian  crossing. 

Description: 

When  the  vehicle  is  in  sight,  the  i  and  m  direction  cosines 
are  tested  and  used  to  interpolate  for  the  times  of  polar  and  meri¬ 
dian  crossings.  When  the  acquisition  time  is  to  be  determined,  the 
integration  routines  CINT  and  EINT  or  the  two-body  routine  KEPLER  is 
employed  to  predict  ahead  until  the  vehicle  comes  into  view. 

Called  By:  OBSERA 

Subroutines  Called;  CINT,  CRSTRE,  EINT,  ERSTRE,  KEPLER 

CCHKEF 

Purpose; 

This  subroutine  tests  criteria  for  changing  the  reference  body 
when  the  Cowell  integrator  la  used.  Even  though  a  change  of  reference 
body  in  the  Cowell  method  is  not  necessary,  it  is  utilised  in  the 
program  in  the  same  manner  as  in  the  Encke  integrator# 
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Descrigtljonj 

Criteria  for  changing  the  reference  body  based  on  regions  of 
Influence  are  tested. 


Called  By:  CITGRA 

Subroutines  Called:  DDOT,  EPHEM,  SERVCE 


CINT  (EINT) 

Purpose : 

This  subroutine  performs  the  integration  when  the  Cowell  formu¬ 
lation  Is  employed. 

Description} 

The  Gill  modification  of  Rungc-Kutta  Is  used  for  short-term  Inte¬ 
gration  and  as  a  starting  procedure  for  the  Nordsleck  long-term  Inte¬ 
gration. 

Called  By:  ATIM,  CITGRA,  MAINA,  OBSERA 
Subroutines  Called  I  CCHREF,  CINT,  CRSTRE 

CITGRA 

Purpose: 

This  subroutine  serves  as  the  sub-main  program  governing  calls 
to  the  Integration  routines  In  the  Cowell  method. 

Description: 

If  not  In  powered  flight,  the  program  checks  for  a  change  of  re¬ 
ference  bodies.  The  delta  of  Integration  Is  determined  by  the  dis¬ 
tance  of  the  vehicle  from  the  refetence  body  and  Integration  Is  per¬ 
formed  up  to  the  next  activity  time  by  calling  the  Integration  routine. 
Called  By:  MAINA 

^ubrou£lnj^^£^jg^:  CCHREF,  CINT,  CRSTRE 

CRSTRE  (ICR) 

Putyose: 

This  subroutine  saves  or  restbres  time,  position,  and  velocity 
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and  acceleration  of  the  vehicle  at  any  designated  time  for  the  Cowell 
integrator. 

Description; 

Depending  on  the  value  of  the  argument  (IRC)  the  information  is 
saved  (IRC  ^  1)  or  restored  (IRC  >  1), 

Called  By;  ATIM,  CITGRA,  MAINA,  OBSERA 

DDOT(A.B) 

Purpose : 

This  function  computes  the  dot  product  of  two  vectors. 
Description; 

The  program  uses  the  formula 

DDOT  -  A(l)*  B(l)  +  A(2)*  B(2)  +  A(3)*  B(3) 

Called  By;  CCHREF,  DERIV,  ECHREF,  LUNOBL,  MVDRAG,  OBD, 

OBLDRG,  OBSERA,  PRINTA,  RECT,  STACUL 

DENSTY(X1,  X2,  X3,  DENP,  DENEQ) 

Purpose; 

This  subroutine  evaluates  the  density  of  air  in  the  upper 
atmosphere. 

Description; 

Tables  of  the  logarithm  of  p  •  C^^  obtained  from  the  Harris- 
Priester  data  are  stored  in  this  routine.  Interpolation  is  performed 
to  calculate  log(p  •  ^av)  at  the  equator  as  a  function  of  the 

altitude  of  the  vehicle  (XI),  solar  flux  (X2) ,  and  solar  time  (X3) • 
Interpolation  is  also  performed  to  calculate  log(p  *  Cav)(DENP)  at 
the  poles,  as  a  function  of  the  altitude  and  solar  flux. 

Called  By;  OBLDRG 

DERIV 

Purpose; 

This  subroutine  evaluates  the  acceleration  terms  for  the  Encke 
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and  Cowell  Integrators* 

Description; 

This  subroutine  computes  the  planetary  perturbations  and  the  solar 
radiation  pressure  perturbations  as  presented  in  the  theory  section  of 
this  manual*  Earth  oblateness  perturbations  and  drag  are  computed  by 
calling  OBLDRG*  Provisions  are  Included  for  computing  powered  flight 


accelerations* 

Called  By;  CINT^  EINT 

Subroutines  Called;  DDOT^  DOMUD^  ECLIPS^  EPHEM^  KEPLER^ 

LUNOBL,  MVDRAG^  OBLDRG^  PFLGHT,  SERVCE 

DMTML  (A^  B>  C,  I,  Jg  K,  lAC,  JAB>  KBC,  IFLAG) 

Purpose; 


This  subroutine  multiplies  two  matrices  of  any  size  (up  to  26  ^ 
26),  or  two  matrices  in  which  the  second  is  transposed* 

Description; 

The  matrices  A,  B,  and  C  are  dimensioned  I  x  J,  K  x  L,  and 
L  X  M  respectively* 


When  IFLAG  *  0,  1,  multiplication  is  done  by  rows  of  A  so 
that  A  can  be  overwritten  if  desired  (l*e*,  the  argument  C  in  the 
calling  sequence  may  actually  be  the  same  as  A) * 

When  IFLAG  *  2,  3,  multiplication  is  done  by  columns  of  B  and 
the  result  is  stored  in  B* 

The  following  table  summarizes  the  results  of  the  program  for 
various  values  of  IFLAG* 


IFLAG 

0 

1 

2 

3 


COMPUTE 

T 

A  •  B 
A  •  B 

T 

A  •  B 
A  •  B 


STORE  RESULT  IN 

C 

C 

B 

B 


lAC  is  the  number  of  rows  of  A  and  C  to  be  used, 

JAB  is  the  number  of  columns  of  A  and  rows  of  B 
(or  B*^)  to  be  used. 
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IBC  is  the  number  of  columns  of  B  (or  bI^)  and 
C  to  be  used* 

Called  By:  LUNOBL,  NUTPRE,  OBLDRG,  OBSERA, 

STACUL,  STAPOS,  XFORM 

DOMUD  (TEST) 

Purpose : 

This  subroutine  checks  for  errors. 

Description; 

The  Overflow  and  Divide  Check  indicators  are  checked.  If  either 

is  on,  and  the  argument  TEST  is  not  equal  to  zero,  TEST  is  assumed  to 

be  a  BCD  word  and  the  message  "ERROR  IN  'TEST'"  is  printed. 

Called  By:  DERIV,  KEPLER,  NUTPRE,  OBLDRG,  OBSERA, 

PRINT A,  RECT,  STAPOS,  XFORM 

ECHREF 

Purpose; 

This  subroutine  tests  criteria  for  changing  the  reference  body 
when  the  Encke  Integrator  is  used. 

Description; 

Criteria  for  changing  the  reference  body  based  on  regions  of 
Influence  are  tested. 

Called  By;  EITGRA 

Subroutines  Called;  DDOT,  EPHEM,  KEPLER,  SERVCE 

ECLIPS  (J,  K) 

Purpose ; 

This  subroutine  checks  for  a  change  in  the  illumination  of  the 
vehicle  and  prints  out  an  appropriate  message  when  the  vehicle 
changes  between  the  three  states:  umbra,  penumbra,  or  sunlight. 
Description; 

This  first  argument  in  the  calling  sequence  indicates  whether 
the  reference  body  is  the  Earth  (J«l) ,  the  Moon  (J*2) ,  or  some  other 
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body  (J>2).  Thtt  second  argument  indicates  whether  the  vehicle  la  In 
the  umbra  (K-1) ,  penumbra  (K»2) ^  or  sunlight  (K*3)#  Values  of  the 
second  argument  are  saved  from  one  entry  to  the  next  and  hence  a 
simple  comparison  is  made  for  determining  a  change  in  vehicle  status* 

If  a  change  has  occurred »  a  message  is  printed  containing  the  refer¬ 
ence  body,  the  time  of  change  and  the  entered  state  (umbra,  penumbra, 
sunlight) • 

Called  By:  DERIV 

EINT  (lENT) 

Purpose; 

This  subroutine  performs  the  integration  when  the  Encke  method 
is  employed* 

Description; 

The  Gill  modification  of  Runge-Kutta  is  used  for  short-term  in¬ 
tegration  and  as  a  starting  procedure  for  the  Nordsieck  long-term 
integration. 

Called  By:  DERIV 

Subroutines  Called;  ATIM,  EITGRA,  OBSERA 

EITGRA 

Purpose; 

This  subroutine  serves  as  the  sub-main  program  governing  calls 
to  the  integration  routines  in  the  Encke  method* 

Description; 

If  not  in  powered  flight,  the  program  checks  for  a  change  of  re¬ 
ference  bodies*  The  delta  of  integration  is  determined  by  the  distance 
of  the  vehicle  from  the  reference  body  and  integration  is  performed  up 
to  the  next  activity  time  by  calling  the  integration  routine*  Fre¬ 
quent  checks  are  also  made  on  the  magnitude  of  the  deviations  from  the 
nominal  orbit*  If  the  magnitudes  increase  beyond  specified  limits,  a 
rectification  of  the  orbit  is  accomplished  by  calling  RECT* 
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Called  By;  MAINA 


Subroutines  Called;  ECHREF,  EIHT,  ERSTRE,  KEPLER,  RECT 

EPHEM 

Purpose : 

This  subroutine  evaluates  the  position  and  velocity  vectors  of 
each  of  the  seven  bodies  with  respect  to  the  reference  body* 
Description! 

Tabular  planetary  positions  are  read  from  an  ephemeris  tape  and 
Everett's  interpolation  formula  is  fitted  to  six  points*  It  is  eval¬ 
uated  to  find  the  position  vector  and  Its  derivative  is  evaluated  to 
find  the  velocity. 

Called  By;  CCKREF,  DERIV,  ECHREF,  STaCUL 
Subroutines  Called:  SERVCE 

ERSTRE  (IRC) 

Purpose  t 

This  subroutine  saves  or  restores  time,  position,  velocity,  and 
acceleration  of  the  vehicle  at  any  designated  time  for  the  Encke  in¬ 
tegrator. 

Description; 

Depending  on  the  value  of  the  argument  (IRC)  the  information  is 
saved  (IRC  ^  1)  or  restored  (IRC  >  1). 

Called  By;  ATIM,  EITGRA,  MAINA,  OBSERA 

FIX  (KTEMP,  ITEMP,  KNAME) 

Purpose; 

This  subroutine  unpacks  a  word  into  five  separate  words* 
Description; 

The  argument  KTEMP  is  assumed  to  contain  a  number  of  the  form 
VVWWXXYYZZ*  The  word  is  unpacked  and  stored  as  follows; 

KNAME  -  W 
ITEMP (1)  -  ZZ 
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ITEMP(2)  -  YY 
ITEMP(3)  -  XX 
ITEMP(4)  -  WW 
Called  By:  OBSERA 

INPUTA 

Purpoee; 

This  subroutine  reeds  all  data  necessary  for  one  run. 

Description; 

The  subroutine  Initializes  necessary  data  and  reads  In  sections 
desired. 

Called  By;  EXECA 
Subroutines  Called;  XFORM 

KEPLER 

Purpose; 

This  subroutine  computes  the  tvo^body  position  and  velocity 
vectors. 

Description; 

A  Newton  Raphson  scheme  Is  used  to  determine  the  differential  ec¬ 
centric  anomaly.  After  convergence,  the  tvorbody  position  end  velo¬ 
city  vectors  are  evalueted.  The  subroutine  uses  Herrick's  method. 
Called  By;  ATIM,  DERIV,  ECHREF,  EITGRA,  OBSERA,  STACUL 
Subroutines  Called;  SERVCE,  DOMUD 

LUNOBL(K) 

Purpose; 

This  subroutine  computes  the  acceleration  due  to  lunar  oblateness. 
Optionally,  It  can  compute  the  llbratlon  and  effect  of  the  eerth's 
J20  term. 

Description; 

When  K  1,  the  llbratlon  matrix  Is  computed  and  then  precessed 
and  nutated. 


ISl 


Whan  K  -  2,  tha  Earth’s  J20  oblatanaaa  tarn  la  calculatadt 
When  K  *  3,  both  lunar  and  Xarth  oblatanaaa  (J20  tarn)  art 
computed. 

Called  Byt  DERIV,  OBD 

Subroutines  Called;  DDOT,  DMIML,  NUTPRE,  SERVCE 

MAINA 

Purpose t 

This  subroutine  handles  the  main  flow  of  the  prograa. 
DeacrJ^^l^^: 

The  routine  calls  subroutines  which  datariBlna  the  next  activity 
time,  govern  the  iijtegration,  compute  tha  observations,  and  print  the 
chosen  Information. 

Called  Byt  EXECA 

Subroutines  Called >  CITGRA,  CRSTRl,  IITCRA,  ERSTRE,  OBSEIA, 

PFINIT,  PRINTA,  lECT,  SERVCE,  STACUL, 

TIMNGA 


MODELA(K) 

Purpose} 

This  subroutine  computes  the  Index  of  refraction  for  tha  tropo¬ 
sphere  or  Ionosphere. 

Description; 

The  index  of  refraction  la  coaputad  with  tha  tropospheric  aodal 
when  K  1  and  the  ionospheric  model  when  K  ■  2. 

Called  By>  OBSERA 


MVP  RAG 


Purpose; 

This  subroutine  computes  the  fsrturbetions  due  to  drag  in  the 
atmospheres  of  Mars  or  Venus. 

Description; 

The  coefficient  of  dreg  is  deteminod  by  interpoletion  froa  given 


152 


tables  as  a  function  of  altitude  and  the  velocity  of  the  vehicle. 

Called  Byt  DERIV 
SubrayJ;lr^eg__Calle^i  DDOT,  SERVCE 

NUTPRE(K) 

Purpose t 

This  subroutine  computes  the  expressions  for  determining  the 
Greenwich  hour  angle,  as  well  as  the  rotational  matrices  between  co¬ 
ordinate  systems. 

Description; 

When  K  1,  the  expressions  used  In  the  nutation  matrix  and  the 
llbratlon  matrix  are  computed. 

When  K  -  2,  the  rotation  matrix  through  the  Greenwich  hour 
angle  Is  computed. 

When  K  3,  the  precession-nutation  matrix  Is  computed. 

Called  Bv!  LUNOBL,  OBLDRG,  PRINTA,  STAPOS,  XFORM 
Subroutines  Celled l  DMTML,  DOMUD 

OBD 

Purpose; 

This  subroutine  computes  observations  from  on-board  Instrtioen- 
tatlon. 

Description; 

The  subroutine  uses  present  position  of  vehicle  with  respect  to 
reference  bodies,  landmarks,  or  ground  stations  to  determine  obser¬ 
vations. 

Called  By;  OBSERA 

Subroutines  Called;  DDOT,  LUNOBL,  SERVCE,  STAPOS 

OBLDRG 

Purpose; 

This  subroutine  computes  the  oblateness  and  air  drag  perturbations 
due  to  the  Earth  and  Its  atmosphere. 
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Description: 


Oblateness  is  computed  in  a  flexible  algorithm  which  can  be  used 
for  zonal,  tesseral,  or  sectorial  harmonics  of  any  order.  The  drag 
perturbation  is  computed  from  the  Harris-Priester  tables  (at  high  alti¬ 
tudes)  and  the  U.  S,  Standard  Atmosphere  (at  low  altitudes). 

Called  By:  DERIV 

Subroutines  Called:  DDOT,  DENSTY,  DKTML,  DOMUD,  NUTPRE,  SERVCE 

OBSERA 

Purpose: 

This  subroutine  computes  the  observables  as  seen  from  a  given 
ground  station. 

Description; 

Given  present  vehicle  location  and  ground  station  location,  both 
referenced  to  inertial  space,  the  observed  values  are  computed.  When 
specified,  refraction  corrections  are  included.  When  in  mode  A,  the 
total  time  in  sight  at  the  station  is  also  computed. 

Called  By:  MAINA 

Subroutines  Called:  ATIM,  CINT,  CRSTRE,  DDOT,  DMTML,  DOMUD, 

EINT,  ERSTRE,  FIX,  KEPLER,  MODELA,  OBD, 

SERVCE,  STAPOS 

PFINIT 

Purpose; 

This  subroutine  performs  initialization  procedures  at  the  entry 
of  a  burn  period. 

Description; 

The  subroutine  determines  coefficients  of  six  Chebyshev  polyno¬ 
mials  which  are  valid  for  the  first  burn  period. 

Called  By;  MAINA 
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PFLCHT 

Purpose: 


This  routine  computes  the  effect  of  a  thrust  perturbation  In  the 
Encke  method. 

Description: 

The  thrust  acceleration  vector  at  the  Initial  burn  time,  the  vehi¬ 
cle  mass  and  mass  rate  are  converted  to  a  set  of  coefficients  for  6 
Chabyshev  polynomials.  These  coefficients  are  determined  In  PFINIT. 

The  subroutine,  PFLGHT,  uses  these  coefficients  to  describe  the  powered 
flight  trajectory  as  a  function  of  time.  Effectively,  the  subroutine 
replaces  KEPLER  In  Encke 's  method.  Integration  of  the  equations  of 
motion  continues  exactly  as  If  powered  flight  was  not  Involved. 

Called  By:  DERIV 
Subroutlnet  Celled ;  SERVCE 

PRINTA 

PurpoM: 

This  subroutine  prints  out  current  trajectory  Information  and 
observations. 

Description: 

Depending  on  whether  the  time  entered  corresponds  to  a  print  ‘time, 
the  program  checks  each  of  the  ten  elements  In  the  IPSEC  array.  If  the 
value  Is  >  0,  the  corresponding  section  Is  printed. 

To  determine  whether  It  Is  a  print  time,  the  subroutine  first 
checks  to  see  whether  the  present  time  Is  within  the  print  portion 

(DTPI)  of  a  total  print  period  (TAU).  If  not,  no  printing  Is  done. 

If  so.  It  next  checks  the  value  of  the  print  Interval  within  DTPI 
(PRATE).  ,  If  It  Is  negative.  It  always  prints.  If  It  Is  positive, 
and  It  Is  the  first  time  Into  the  present  print  period.  It  prints, 

otherwise  no  printing  Is  done.  When  MDE  -  1  or  4,  printing  occurs 

at  each  entry  to  the  routine. 
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Called  By:  MAINA 

Subroutlneg  Called;  DDOT,  DOMUD,  NUTPRE,  SERVCE 


RECT 

Purpose; 

This  subroutine  computes  the  parameters  pertinent  to  a  rectifi¬ 
cation  in  Encke's  method. 

Description; 

The  two-body  position  and  velocity  vectors  are  equated  to  the 
true  position  and  veolcity  vectors.  In  addition,  these  components  are 
saved.  Perturbations  in  position  and  velocity  are  set  equal  to  zero, 
and  elements  used  by  the  KEPLER  subroutine  are  computed. 

Called  By;  EITGRA,  MAINA 

Subroutines  Called;  DDOT ,  DOMUD 

SERVCE  (A,  B,  C,  I) 

Purpose; 

This  subroutine  computes  the  cross  product  of  two  vectors  as  well 
as  the  magnitude,  magnitude  squared  and  magnitude  cubed  of  a  vector. 

Description: 

When  1*1,  the  cross  product  of  A  and  B  is  computed  and 
the  result  stored  in  C(l),  C(2),  and  C(3),  and  then  continues  as 
when  I  •  2, 

When  I  ■  2,  the  magnitude  cubed,  magnitude,  and  magnitude 
squared  is  stored  in  C(4),  C(5),  and  C(6),  respectively. 

Called  By:  CCHREF,  DERIV,  ECHREF,  EPHEM,  KEPLER,  LUNOBL, 

MAINA,  MVDRAG,  OBD,  OBLDRG,  OBSERA,  PFLGHT, 

PRINTA,  STACUL 


STACUL 

Purpose: 

This  subroutine  determines  the  time  of  occultatlon  (YCOM(23))« 
One  of  two  types  of  occultatlon  may  be  considered,  vehicle  or  star 
occultatlon* 
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Description: 

A  Newton-Raphson  Iteration  la  used  on  the  two-body  solution  to 
determine  the  time  of  occultatlon.  When  considering  star  occultatlon, 
the  occulting  body  Is  the  reference  body,  except  In  the  Earth-Moon 
system.  In  which  either  the  Moon  or  the  Earth  can  do  the  occulting. 
Vehicle  occultatlon  occurs  only  In  Moon  reference.  Up  to  three  dif¬ 
ferent  ground  stations  may  be  considered  for  occultatlon. 

Called  By:  MAINA 

Subroutines  Called:  DDOT,  DMTML,  EPHEM,  KEPLER,  SERVCE,  STAPOS 

STAPOS 

Purpose: 

This  subroutine  computes  the  station  position  and  velocity  vectors. 

Description: 

Given  the  geodetic  latitude,  longitude,  and  altitude  of  the 
station,  the  coordinates  of  the  station  In  the  Greenwich  coordinate 
system  are  computed.  A  rotation  through  the  Greenwich  hour  angle  then 
brings  the  coordinates  Into  the  Inertial  system.  The  velocity  of  the 
station  Is  computed  from  the  rotational  speed  of  the  earth. 

Called  By:  OBD,  OBSERA,  STACUL 

Subroutines  Called:  DMTML,  DOMUD,  NUTPRE 

TIMNGA 

Purpose: 

This  subroutine  determines  the  next  activity  time  (l.e.,  the  next 
time  at  which  an  observation  Is  to  be  computed  or  a  printout  of  the 
trajectory  Is  to  occur). 

Description: 

Logic  Is  set  up  for  establishing  an  array  of  times  of  Interest 
from  which  the  earliest  time  Is  selected.  Flags  are  set  when  the 
time  selected  Is  the  final  time  or  when  e  time  Is  repeated. 
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Callad  By;  MAINA 


XFORM 

Purpose; 

This  routine  accepts  the  Input  Information  concerning  the  vehi¬ 
cle's  Initial  conditions,  and  transforms  them  Into  the  proper  units 
and  coordinate  system. 

Description; 

After  the  Initial  conditions  and  associated  flags  are  read  In 
(Section  2  of  the  group  I  Inputs),  the  program  converts  them  into 
the  position  and  velocity  vectors  In  the  units  of  the  Earth  radii, 
and  Earth  radll/hour.  Then,  If  desired,  the  vectors  are  transformed 
Into  the  base  date  system. 

Called  By;  INPUTA 

Subroutines  Called;  DMTML,  DOMUD,  NUTPRE 
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